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.

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

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

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

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

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

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

          INFO is INTEGER
          = 0:  successful exit

 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.

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

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

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

 ILAZLC scans A for its last non-zero column.

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

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

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

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

 ILAZLR scans A for its last non-zero row.

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

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

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

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

 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.

          N is INTEGER
          The number of components of the vector x.

          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 is COMPLEX*16 array, dimension
                         (1+(N-1)*abs(INCX))
          The n-element vector x.

          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

 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

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

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

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

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).

          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 is DOUBLE PRECISION array, dimension (NB)
          The off-diagonal elements of the first NB rows and columns of
          the reduced matrix.

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

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

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

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

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

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

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

 ZLACGV conjugates a complex vector of length N.

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

          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 is INTEGER
          The spacing between successive elements of X.

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

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

          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 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 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 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 is INTEGER array, dimension (3)
         ISAVE is used to save variables between calls to ZLACN2

  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)

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

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

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

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

          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 is INTEGER
          The number of rows of the matrix A.  M >= 0.

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

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).

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

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

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

          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 is INTEGER
          The number of rows of the matrix A.  M >= 0.

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

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).

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

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

 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.

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

          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 is COMPLEX*16 array, dimension (LDA, N)
          On entry, A contains the M by N matrix A.

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

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

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

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

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

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

 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.

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

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

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

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

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

          C is COMPLEX*16

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

 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.

          X is COMPLEX*16

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

 ZLAEIN uses inverse iteration to find a right or left eigenvector
 corresponding to the eigenvalue W of a complex upper Hessenberg
 matrix H.

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

          NOINIT is LOGICAL
          = .TRUE. : no initial vector supplied in V
          = .FALSE.: initial vector supplied in V.

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

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

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

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

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

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

          RWORK is DOUBLE PRECISION array, dimension (N)

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

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

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

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

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

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

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

          RT1 is DOUBLE PRECISION
         The eigenvalue of larger absolute value.

          RT2 is DOUBLE PRECISION
         The eigenvalue of smaller absolute value.

          CS1 is DOUBLE PRECISION

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

  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.

 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.

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

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

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

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

          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 is INTEGER
          The leading dimension of the array SA.  LDSA >= max(1,M).

          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.

 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.

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

          A1 is DOUBLE PRECISION

          A2 is COMPLEX*16

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

          B1 is DOUBLE PRECISION

          B2 is COMPLEX*16

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

          CSU is DOUBLE PRECISION

          SNU is COMPLEX*16
          The desired unitary matrix U.

          CSV is DOUBLE PRECISION

          SNV is COMPLEX*16
          The desired unitary matrix V.

          CSQ is DOUBLE PRECISION

          SNQ is COMPLEX*16
          The desired unitary matrix Q.

 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.

          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 is INTEGER
          The order of the matrix A.  N >= 0.

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

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

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

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

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

          X is COMPLEX*16 array, dimension (LDX,NRHS)
          The N by NRHS matrix X.

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

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

          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 is INTEGER
          The leading dimension of the array B.  LDB >= max(N,1).

    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.

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

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

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

          ILO is INTEGER

          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 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 is INTEGER
          The leading dimension of the array H. LDH >= max(1,N).

          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 is INTEGER

          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 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 is INTEGER
          The leading dimension of the array Z. LDZ >= max(1,N).

          INFO is INTEGER
           =   0: successful exit
          .GT. 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 .GT. 0 and WANTT is .FALSE., then on exit,
                  the remaining unconverged eigenvalues are the
                  eigenvalues of the upper Hessenberg matrix
                  rows and columns ILO thorugh INFO of the final,
                  output value of H.

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

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

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

 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.

          N is INTEGER
          The order of the matrix A.

          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 is INTEGER
          The number of columns to be reduced.

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).

          TAU is COMPLEX*16 array, dimension (NB)
          The scalar factors of the elementary reflectors. See Further
          Details.

          T is COMPLEX*16 array, dimension (LDT,NB)
          The upper triangular matrix T.

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

          Y is COMPLEX*16 array, dimension (LDY,NB)
          The n-by-nb matrix Y.

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

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

 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.

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

          J is INTEGER
          Length of X and W

          X is COMPLEX*16 array, dimension (J)
          The j-vector x.

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

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

          GAMMA is COMPLEX*16
          The diagonal element gamma.

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

          S is COMPLEX*16
          Sine needed in forming xhat.

          C is COMPLEX*16
          Cosine needed in forming xhat.

 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.

    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.

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

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

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

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

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

 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.

    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.

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

          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 is INTEGER
          The order of the matrix A.  N >= 0.  When N = 0, ZLANHB is
          set to zero.

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

          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 is INTEGER
          The leading dimension of the array AB.  LDAB >= K+1.

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

 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.

    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.

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

          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 is INTEGER
          The order of the matrix A.  N >= 0.  When N = 0, ZLANHP is
          set to zero.

          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 is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
          WORK is not referenced.

 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.

    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.

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

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

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(N,1).

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

 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.

    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.

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

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

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

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

 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.

    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.

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

          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 is INTEGER
          The order of the matrix A.  N >= 0.  When N = 0, ZLANSB is
          set to zero.

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

          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 is INTEGER
          The leading dimension of the array AB.  LDAB >= K+1.

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

 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.

    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.

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

          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 is INTEGER
          The order of the matrix A.  N >= 0.  When N = 0, ZLANSP is
          set to zero.

          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 is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
          WORK is not referenced.

 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.

    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.

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

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

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

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

          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 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 is INTEGER
          The leading dimension of the array AB.  LDAB >= K+1.

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

 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.

    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.

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

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

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

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

          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 is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
          where LWORK >= N when NORM = 'I'; otherwise, WORK is not
          referenced.

 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.

    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.

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

          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 is CHARACTER*1
          Specifies whether or not the matrix A has unit diagonal.
          = 'N':  Non-unit diagonal
          = 'U':  Unit diagonal

          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 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 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 is INTEGER
          The leading dimension of the array A.  LDA >= max(M,1).

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

 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.

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

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

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

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

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

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

 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.

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

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

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

          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 is INTEGER
          The leading dimension of the array X, LDX >= MAX(1,M).

          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.

 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.

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

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

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

          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 is INTEGER
          The leading dimension of the array X, LDX >= MAX(1,M).

          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.

 ZLAQHB equilibrates a Hermitian band matrix A
 using the scaling factors in the vector S.

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 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 is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.

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

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

          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix entry.

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

  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.

 ZLAQHP equilibrates a Hermitian matrix A using the scaling factors
 in the vector S.

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 is DOUBLE PRECISION array, dimension (N)
          The scale factors for A.

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

          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix entry.

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

  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.

 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.

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

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

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

          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 is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).

          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 is COMPLEX*16 array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.

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

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

          WORK is COMPLEX*16 array, dimension (N)

 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.

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

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

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

          NB is INTEGER
          The number of columns to factorize.

          KB is INTEGER
          The number of columns actually factorized.

          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 is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).

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

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

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

          AUXV is COMPLEX*16 array, dimension (NB)
          Auxiliar vector.

          F is COMPLEX*16 array, dimension (LDF,NB)
          Matrix F**H = L * Y**H * A.

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

    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.

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

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

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

          ILO is INTEGER

          IHI is INTEGER

           It is assumed that H is already upper triangular in rows
           and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
           previous call to 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.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
           If N = 0, then ILO = 1 and IHI = 0.

          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.GT.0 is given under the
           description of INFO below.)

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

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

          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 is INTEGER

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

          Z 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.GT.0 is given under
           the description of INFO below.)

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

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

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

           If LWORK = -1, then 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 is INTEGER
             =  0:  successful exit
           .GT. 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 .GT. 0 and WANT is .FALSE., then on exit,
                the remaining unconverged eigenvalues are the eigen-
                values of the upper Hessenberg matrix rows and
                columns ILO through INFO of the final, output
                value of H.

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

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

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

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

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

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

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

      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.

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

          H is COMPLEX*16 array of dimension (LDH,N)
              The 2-by-2 or 3-by-3 matrix H in (*).

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

          S1 is COMPLEX*16

          S2 is COMPLEX*16

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

          V is COMPLEX*16 array of dimension N
              A scalar multiple of the first column of the
              matrix K in (*).

    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.

          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 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 is INTEGER
          The order of the matrix H and (if WANTZ is .TRUE.) the
          order of the unitary matrix Z.

          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 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 is INTEGER
          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).

          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 is integer
          Leading dimension of H just as declared in the calling
          subroutine.  N .LE. LDH

          ILOZ is INTEGER

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

          Z 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 is integer
          The leading dimension of Z just as declared in the
          calling subroutine.  1 .LE. LDZ.

          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 is integer
          The number of converged eigenvalues uncovered by this
          subroutine.

          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 is COMPLEX*16 array, dimension (LDV,NW)
          An NW-by-NW work array.

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

          NH is integer scalar
          The number of columns of T.  NH.GE.NW.

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

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

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

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

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

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

    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.

          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 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 is INTEGER
          The order of the matrix H and (if WANTZ is .TRUE.) the
          order of the unitary matrix Z.

          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 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 is INTEGER
          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).

          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 is integer
          Leading dimension of H just as declared in the calling
          subroutine.  N .LE. LDH

          ILOZ is INTEGER

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

          Z 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 is integer
          The leading dimension of Z just as declared in the
          calling subroutine.  1 .LE. LDZ.

          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 is integer
          The number of converged eigenvalues uncovered by this
          subroutine.

          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 is COMPLEX*16 array, dimension (LDV,NW)
          An NW-by-NW work array.

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

          NH is integer scalar
          The number of columns of T.  NH.GE.NW.

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

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

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

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

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

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

    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.

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

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

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

          ILO is INTEGER

          IHI is INTEGER
           It is assumed that H is already upper triangular in rows
           and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
           previous call to 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.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
           If N = 0, then ILO = 1 and IHI = 0.

          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.GT.0 is given under the
           description of INFO below.)

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

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

          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 is INTEGER

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

          Z 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.GT.0 is given under
           the description of INFO below.)

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

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

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

           If LWORK = -1, then 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 is INTEGER
             =  0:  successful exit
           .GT. 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 .GT. 0 and WANT is .FALSE., then on exit,
                the remaining unconverged eigenvalues are the eigen-
                values of the upper Hessenberg matrix rows and
                columns ILO through INFO of the final, output
                value of H.

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

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

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

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

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

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

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

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

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

          WANTZ is logical scalar
             WANTZ = .true. if the unitary Schur factor is being
             computed.  WANTZ is set to .false. otherwise.

          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: ZLAQR5 accumulates reflections, uses matrix-matrix
             multiply to update the far-from-diagonal matrix entries,
             and takes advantage of 2-by-2 block structure during
             matrix multiplies.

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

          KTOP is integer scalar

          KBOT is integer scalar
             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 is integer scalar
             NSHFTS gives the number of simultaneous shifts.  NSHFTS
             must be positive and even.

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

          H is COMPLEX*16 array of size (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 is integer scalar
             LDH is the leading dimension of H just as declared in the
             calling procedure.  LDH.GE.MAX(1,N).

          ILOZ is INTEGER

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

          Z is COMPLEX*16 array of size (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 is integer scalar
             LDA is the leading dimension of Z just as declared in
             the calling procedure. LDZ.GE.N.

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

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

          U is COMPLEX*16 array of size
             (LDU,3*NSHFTS-3)

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

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

          WH is COMPLEX*16 array of size (LDWH,NH)

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

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

          WV is COMPLEX*16 array of size
             (LDWV,3*NSHFTS-3)

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

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

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 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 is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.

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

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

          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix entry.

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

  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.

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

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 is DOUBLE PRECISION array, dimension (N)
          The scale factors for A.

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

          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix entry.

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

  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.

 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.

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

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

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

          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 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 is DOUBLE PRECISION array, dimension (N)
           The n diagonal elements of the diagonal matrix D.

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

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

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

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

          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 is LOGICAL
           Specifies whether NEGCNT has to be computed.

          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 is DOUBLE PRECISION
           The square of the 2-norm of Z.

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

          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 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 is DOUBLE PRECISION
           NRMINV = 1/SQRT( ZTZ )

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

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

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

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

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

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

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

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

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

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

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

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

 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.

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

          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 is DOUBLE PRECISION array, dimension (LDA, M)
          On entry, A contains the M by M matrix A.

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

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

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

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

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

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

 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.

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

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

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

          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 is INTEGER
          The increment between elements of v. INCV <> 0.

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

          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 is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is COMPLEX*16 array, dimension
                         (N) if SIDE = 'L'
                      or (M) if SIDE = 'R'

 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.

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

          TRANS is CHARACTER*1
          = 'N': apply H (No transpose)
          = 'C': apply H**H (Conjugate transpose)

          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 is CHARACTER*1
          Indicates how the vectors which define the elementary
          reflectors are stored:
          = 'C': Columnwise
          = 'R': Rowwise

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

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

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

          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 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 is COMPLEX*16 array, dimension (LDT,K)
          The triangular K-by-K matrix T in the representation of the
          block reflector.

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

          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 is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

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

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

  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 )

 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 .

          N is INTEGER
          The order of the elementary reflector.

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

          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 is INTEGER
          The increment between elements of X. INCX > 0.

          TAU is COMPLEX*16
          The value tau.

 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.

          N is INTEGER
          The order of the elementary reflector.

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

          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 is INTEGER
          The increment between elements of X. INCX > 0.

          TAU is COMPLEX*16
          The value tau.

 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

          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 is CHARACTER*1
          Specifies how the vectors which define the elementary
          reflectors are stored (see also Further Details):
          = 'C': columnwise
          = 'R': rowwise

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

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

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

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

          TAU is COMPLEX*16 array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i).

          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 is INTEGER
          The leading dimension of the array T. LDT >= K.

  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 )

 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.

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

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

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

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

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

          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 is INTEGER
          The leading dimension of the array C. LDA >= max(1,M).

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

 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.

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

          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 is INTEGER
          The increment between elements of X. INCX > 0.

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

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

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

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

  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.

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

          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 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 is INTEGER
          The number of random numbers to be generated.

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

  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.

 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.

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

          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 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 is DOUBLE PRECISION array, dimension (N)
          On entry, the N diagonal elements of the diagonal matrix D.
          On exit, D may be overwritten.

          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 is DOUBLE PRECISION
          The minimum pivot allowed in the Sturm sequence.

          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 is INTEGER
          The total number of input eigenvalues.  0 <= M <= N.

          DOL is INTEGER

          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 is DOUBLE PRECISION

          RTOL1 is DOUBLE PRECISION

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

          W 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 is DOUBLE PRECISION array, dimension (N)
          The first M elements contain the semiwidth of the uncertainty
          interval of the corresponding eigenvalue in W

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

          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 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 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 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 is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).

          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 is DOUBLE PRECISION array, dimension (12*N)

          IWORK is INTEGER array, dimension (7*N)

          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.

 ZLARTG generates a plane rotation so that

    [  CS  SN  ]     [ F ]     [ R ]
    [  __      ]  .  [   ]  =  [   ]   where CS**2 + |SN|**2 = 1.
    [ -SN  CS  ]     [ G ]     [ 0 ]

 This is a faster version of the BLAS1 routine ZROTG, except for
 the following differences:
    F and G are unchanged on return.
    If G=0, then CS=1 and SN=0.
    If F=0, then CS=0 and SN is chosen so that R is real.

          F is COMPLEX*16
          The first component of vector to be rotated.

          G is COMPLEX*16
          The second component of vector to be rotated.

          CS is DOUBLE PRECISION
          The cosine of the rotation.

          SN is COMPLEX*16
          The sine of the rotation.

          R is COMPLEX*16
          The nonzero component of the rotated vector.

  3-5-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.

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

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

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

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

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

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

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

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

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

 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.

          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 is INTEGER
          The lower bandwidth of A.  Referenced only if TYPE = 'B',
          'Q' or 'Z'.

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

          CFROM is DOUBLE PRECISION

          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 is INTEGER
          The number of rows of the matrix A.  M >= 0.

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

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

          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 is INTEGER
          0  - successful exit
          <0 - if INFO = -i, the i-th argument had an illegal value.

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

          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 is INTEGER
          On entry, M specifies the number of rows of A.

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

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

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

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).

 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.

          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 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 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 is INTEGER
          The number of rows of the matrix A.  If m <= 1, an immediate
          return is effected.

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

          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 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 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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).

 ZLASSQ returns the values scl and ssq such that

    ( scl**2 )*ssq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,

 where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is
 assumed to be at least unity and the value of ssq will then satisfy

    1.0 .le. ssq .le. ( sumsq + 2*n ).

 scale is assumed to be non-negative and scl returns the value

    scl = max( scale, abs( real( x( i ) ) ), abs( aimag( x( i ) ) ) ),
           i

 scale and sumsq must be supplied in SCALE and SUMSQ respectively.
 SCALE and SUMSQ are overwritten by scl and ssq respectively.

 The routine makes only one pass through the vector X.

          N is INTEGER
          The number of elements to be used from the vector X.

          X is COMPLEX*16 array, dimension (N)
          The vector x as described above.
             x( i )  = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.

          INCX is INTEGER
          The increment between successive values of the vector X.
          INCX > 0.

          SCALE is DOUBLE PRECISION
          On entry, the value  scale  in the equation above.
          On exit, SCALE is overwritten with the value  scl .

          SUMSQ is DOUBLE PRECISION
          On entry, the value  sumsq  in the equation above.
          On exit, SUMSQ is overwritten with the value  ssq .

 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.

          N is INTEGER
          The number of columns of the matrix 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 is INTEGER
          The leading dimension of the array A.

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

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

          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)*INCX of IPIV are accessed.
          IPIV(K) = L implies rows K and L are to be interchanged.

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

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

 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.

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

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

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

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

          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 is INTEGER
          The leading dimension of the array SA.  LDSA >= max(1,M).

          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.

 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.

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

          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 is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit triangular

          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 is INTEGER
          The order of the matrix A.  N >= 0.

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

          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 is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.

          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 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 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 is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -k, the k-th argument had an illegal value

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

 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.

          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 is INTEGER
          The number of columns of the matrix 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 is INTEGER
          The leading dimension of the array Z.  LDA >= max(1, N).

          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 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 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 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 is INTEGER array, dimension (N).
          The pivot indices; for 1 <= j <= N, column j of the
          matrix has been interchanged with column JPIV(j).

1995.

 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.

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

          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 is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit triangular

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 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 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 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 is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -k, the k-th argument had an illegal value

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

 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.

          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 is INTEGER
          The order of the matrix A.

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

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).

          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 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 is COMPLEX*16 array, dimension (LDW,NB)
          The n-by-nb matrix W required to update the unreduced part
          of A.

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

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

 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.

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

          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 is CHARACTER*1
          Specifies whether or not the matrix A is unit triangular.
          = 'N':  Non-unit triangular
          = 'U':  Unit triangular

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max (1,N).

          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 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 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 is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -k, the k-th argument had an illegal value

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

 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.

          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 is INTEGER
          The order of the triangular factor U or L.  N >= 0.

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).

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

 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.

          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 is INTEGER
          The order of the triangular factor U or L.  N >= 0.

          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 is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).

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

 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.

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

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

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

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

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

          C is DOUBLE PRECISION

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

 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.

          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 is INTEGER
           On entry, N specifies the order of the matrix A.
           N must be at least zero.
           Unchanged on exit.

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

          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 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 is INTEGER
           On entry, INCX specifies the increment for the elements of
           X. INCX must not be zero.
           Unchanged on exit.

          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 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 is INTEGER
           On entry, INCY specifies the increment for the elements of
           Y. INCY must not be zero.
           Unchanged on exit.

 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.

          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 is INTEGER
           On entry, N specifies the order of the matrix A.
           N must be at least zero.
           Unchanged on exit.

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

          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 is INTEGER
           On entry, INCX specifies the increment for the elements of
           X. INCX must not be zero.
           Unchanged on exit.

          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.

 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.

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

          TRANS is CHARACTER*1
          = 'N': apply H (No transpose)
          = 'C': apply H**H (Conjugate transpose)

          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 is CHARACTER*1
          Indicates how the vectors which define the elementary
          reflectors are stored:
          = 'C': Columns
          = 'R': Rows

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

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

          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 is INTEGER
          The order of the trapezoidal part of V.
          K >= L >= 0.  See Further Details.

          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 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 is COMPLEX*16 array, dimension (LDT,K)
          The triangular K-by-K matrix T in the representation of the
          block reflector.

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

          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 is INTEGER
          The leading dimension of the array A.
          If SIDE = 'L', LDC >= max(1,K);
          If SIDE = 'R', LDC >= max(1,M).

          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 is INTEGER
          The leading dimension of the array B.
          LDB >= max(1,M).

          WORK is COMPLEX*16 array, dimension
          (LDWORK,N) if SIDE = 'L',
          (LDWORK,K) if SIDE = 'R'.

          LDWORK is INTEGER
          The leading dimension of the array WORK.
          If SIDE = 'L', LDWORK >= K;
          if SIDE = 'R', LDWORK >= M.

  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.