CLABRD 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 CGEBRD

M

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

N

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

NB

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

A

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

LDA

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

D

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

E

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

TAUQ

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

TAUP

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

X

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

LDX

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

Y

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

LDY

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLACGV conjugates a complex vector of length N.

N

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

X

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

INCX

INCX is INTEGER
The spacing between successive elements of X.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

N

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

V

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

X

 X is COMPLEX 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 CLACN2 must be
re-called with all the other parameters unchanged.

EST

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

KASE

 KASE is INTEGER
On the initial call to CLACN2, 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 CLACN2, KASE will again be 0.

ISAVE

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Originally named CONEST, dated March 16, 1988.
Last modified:  April, 1999
This is a thread safe version of CLACON, which uses the array ISAVE
in place of a SAVE statement, as follows:
   CLACON     CLACN2
    JUMP     ISAVE(1)
    J        ISAVE(2)
    ITER     ISAVE(3)

Nick Higham, University of Manchester

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

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

N

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

V

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

X

 X is COMPLEX 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 CLACON must be
re-called with all the other parameters unchanged.

EST

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

KASE

 KASE is INTEGER
On the initial call to CLACON, 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 CLACON, KASE will again be 0.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

Nick Higham, University of Manchester

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

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

UPLO

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

UPLO

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

M

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

N

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

A

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

LDA

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

B

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

LDB

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLACRM 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

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

N

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

A

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

LDA

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

B

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

LDB

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

C

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

LDC

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

RWORK

RWORK is REAL array, dimension (2*M*N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLACRT 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

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

CX

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

INCX

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

CY

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

INCY

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

C

C is COMPLEX

S

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLADIV := 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

X is COMPLEX

Y

Y is COMPLEX
The complex scalars X and Y.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

RIGHTV

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

NOINIT

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

N

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

H

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

LDH

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

W

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

V

V is COMPLEX array, dimension (N)
On entry, if NOINIT = .FALSE., V must contain a starting
vector for inverse iteration; otherwise V need not be set.
On exit, V contains the computed eigenvector, normalized so
that the component of largest magnitude has magnitude 1; here
the magnitude of a complex number (x,y) is taken to be
|x| + |y|.

B

B is COMPLEX array, dimension (LDB,N)

LDB

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

RWORK

RWORK is REAL array, dimension (N)

EPS3

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

SMLNUM

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAEV2 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

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

B

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

C

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

RT1

 RT1 is REAL
The eigenvalue of larger absolute value.

RT2

 RT2 is REAL
The eigenvalue of smaller absolute value.

CS1

CS1 is REAL

SN1

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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.

CLAGS2 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

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

A1

A1 is REAL

A2

A2 is COMPLEX

A3

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

B1

B1 is REAL

B2

B2 is COMPLEX

B3

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

CSU

CSU is REAL

SNU

SNU is COMPLEX
The desired unitary matrix U.

CSV

CSV is REAL

SNV

SNV is COMPLEX
The desired unitary matrix V.

CSQ

CSQ is REAL

SNQ

SNQ is COMPLEX
The desired unitary matrix Q.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAGTM 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

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

N

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

NRHS

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

ALPHA

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

DL

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

D

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

DU

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

X

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

LDX

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

BETA

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

B

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

LDB

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAHQR 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

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

WANTZ

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

N

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

ILO

ILO is INTEGER

IHI

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

H

H is COMPLEX array, dimension (LDH,N)
On entry, the upper Hessenberg matrix H.
On exit, if INFO is zero and if WANTT is .TRUE., then H
is upper triangular in rows and columns ILO:IHI.  If INFO
is zero and if WANTT is .FALSE., then the contents of H
are unspecified on exit.  The output state of H in case
INF is positive is below under the description of INFO.

LDH

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

W

W is COMPLEX array, dimension (N)
The computed eigenvalues ILO to IHI are stored in the
corresponding elements of W. If WANTT is .TRUE., the
eigenvalues are stored in the same order as on the diagonal
of the Schur form returned in H, with W(i) = H(i,i).

ILOZ

ILOZ is INTEGER

IHIZ

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

Z

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

LDZ

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

INFO

INFO is INTEGER
 = 0:  successful exit
 > 0:  if INFO = i, CLAHQR failed to compute all the
        eigenvalues ILO to IHI in a total of 30 iterations
        per eigenvalue; elements i+1:ihi of W contain
        those eigenvalues which have been successfully
        computed.
        If INFO > 0 and WANTT is .FALSE., then on exit,
        the remaining unconverged eigenvalues are the
        eigenvalues of the upper Hessenberg matrix
        rows and columns ILO through INFO of the final,
        output value of H.
        If INFO > 0 and WANTT is .TRUE., then on exit
(*)       (initial value of H)*U  = U*(final value of H)
        where U is an orthogonal matrix.    The final
        value of H is upper Hessenberg and triangular in
        rows and columns INFO+1 through IHI.
        If INFO > 0 and WANTZ is .TRUE., then on exit
            (final value of Z)  = (initial value of Z)*U
        where U is the orthogonal matrix in (*)
        (regardless of the value of WANTT.)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

N

N is INTEGER
The order of the matrix A.

K

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

NB

NB is INTEGER
The number of columns to be reduced.

A

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

LDA

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

TAU

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

T

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

LDT

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

Y

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

LDY

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

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

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

J

J is INTEGER
Length of X and W

X

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

SEST

SEST is REAL
Estimated singular value of j by j matrix L

W

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

GAMMA

GAMMA is COMPLEX
The diagonal element gamma.

SESTPR

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

S

S is COMPLEX
Sine needed in forming xhat.

C

C is COMPLEX
Cosine needed in forming xhat.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANGT

   CLANGT = ( 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

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

N

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

DL

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

D

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

DU

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANHB

   CLANHB = ( 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

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

UPLO

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

N

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

K

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

AB

AB is COMPLEX array, dimension (LDAB,N)
The upper or lower triangle of the hermitian band matrix A,
stored in the first K+1 rows of AB.  The j-th column of A is
stored in the j-th column of the array AB as follows:
if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)   = A(i,j) for j<=i<=min(n,j+k).
Note that the imaginary parts of the diagonal elements need
not be set and are assumed to be zero.

LDAB

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANHP

   CLANHP = ( 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

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

UPLO

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

N

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

AP

AP is COMPLEX array, dimension (N*(N+1)/2)
The upper or lower triangle of the hermitian matrix A, packed
columnwise in a linear array.  The j-th column of A is stored
in the array AP as follows:
if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
Note that the  imaginary parts of the diagonal elements need
not be set and are assumed to be zero.

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANHS

   CLANHS = ( 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

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

N

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

A

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

LDA

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANHT

   CLANHT = ( 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

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

N

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

D

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

E

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANSB

   CLANSB = ( 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

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

UPLO

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

N

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

K

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

AB

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

LDAB

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANSP

   CLANSP = ( 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

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

UPLO

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

N

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

AP

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANTB

   CLANTB = ( 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

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

UPLO

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

DIAG

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

N

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

K

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

AB

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

LDAB

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANTP

   CLANTP = ( 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

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

UPLO

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

DIAG

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

N

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

AP

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLANTR

   CLANTR = ( 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

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

UPLO

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

DIAG

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

M

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

N

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

A

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

LDA

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

X

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

INCX

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

Y

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

INCY

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

SSMIN

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAPMR 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

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

M

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

N

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

X

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

LDX

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

K

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAPMT 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

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

M

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

N

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

X

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

LDX

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

K

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAQHB equilibrates an Hermitian band matrix A using the scaling
factors in the vector S.

UPLO

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

N

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

KD

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

AB

AB is COMPLEX array, dimension (LDAB,N)
On entry, the upper or lower triangle of the symmetric band
matrix A, stored in the first KD+1 rows of the array.  The
j-th column of A is stored in the j-th column of the array AB
as follows:
if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
On exit, if INFO = 0, the triangular factor U or L from the
Cholesky factorization A = U**H *U or A = L*L**H of the band
matrix A, in the same storage format as A.

LDAB

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

S

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

SCOND

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

AMAX

AMAX is REAL
Absolute value of largest matrix entry.

EQUED

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

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

UPLO

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

N

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

AP

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

S

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

SCOND

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

AMAX

AMAX is REAL
Absolute value of largest matrix entry.

EQUED

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

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAQP2 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

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

N

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

OFFSET

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

A

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

LDA

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

JPVT

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

TAU

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

VN1

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

VN2

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

WORK

WORK is COMPLEX array, dimension (N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

LAPACK Working Note 176

CLAQPS 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

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

N

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

OFFSET

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

NB

NB is INTEGER
The number of columns to factorize.

KB

KB is INTEGER
The number of columns actually factorized.

A

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

LDA

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

JPVT

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

TAU

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

VN1

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

VN2

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

AUXV

AUXV is COMPLEX array, dimension (NB)
Auxiliary vector.

F

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

LDF

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA

LAPACK Working Note 176

CLAQR0 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

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

WANTZ

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

N

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

ILO

ILO is INTEGER

IHI

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

H

H is COMPLEX array, dimension (LDH,N)
 On entry, the upper Hessenberg matrix H.
 On exit, if INFO = 0 and WANTT is .TRUE., then H
 contains the upper triangular matrix T from the Schur
 decomposition (the Schur form). If INFO = 0 and WANT is
 .FALSE., then the contents of H are unspecified on exit.
 (The output value of H when INFO > 0 is given under the
 description of INFO below.)
 This subroutine may explicitly set H(i,j) = 0 for i > j and
 j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.

LDH

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

W

W is COMPLEX array, dimension (N)
 The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
 in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
 stored in the same order as on the diagonal of the Schur
 form returned in H, with W(i) = H(i,i).

ILOZ

ILOZ is INTEGER

IHIZ

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

Z

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

LDZ

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

WORK

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

LWORK

LWORK is INTEGER
 The dimension of the array WORK.  LWORK >= max(1,N)
 is sufficient, but LWORK typically as large as 6*N may
 be required for optimal performance.  A workspace query
 to determine the optimal workspace size is recommended.
 If LWORK = -1, then CLAQR0 does a workspace query.
 In this case, CLAQR0 checks the input parameters and
 estimates the optimal workspace size for the given
 values of N, ILO and IHI.  The estimate is returned
 in WORK(1).  No error message related to LWORK is
 issued by XERBLA.  Neither H nor Z are accessed.

INFO

INFO is INTEGER
   = 0:  successful exit
   > 0:  if INFO = i, CLAQR0 failed to compute all of
      the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
      and WI contain those eigenvalues which have been
      successfully computed.  (Failures are rare.)
      If INFO > 0 and WANT is .FALSE., then on exit,
      the remaining unconverged eigenvalues are the eigen-
      values of the upper Hessenberg matrix rows and
      columns ILO through INFO of the final, output
      value of H.
      If INFO > 0 and WANTT is .TRUE., then on exit
 (*)  (initial value of H)*U  = U*(final value of H)
      where U is a unitary matrix.  The final
      value of  H is upper Hessenberg and triangular in
      rows and columns INFO+1 through IHI.
      If INFO > 0 and WANTZ is .TRUE., then on exit
        (final value of Z(ILO:IHI,ILOZ:IHIZ)
         =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
      where U is the unitary matrix in (*) (regard-
      less of the value of WANTT.)
      If INFO > 0 and WANTZ is .FALSE., then Z is not
      accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part II: Aggressive Early Deflation, SIAM Journal of Matrix Analysis, volume 23, pages 948--973, 2002.

Given a 2-by-2 or 3-by-3 matrix H, CLAQR1 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

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

H

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

LDH

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

S1

S1 is COMPLEX

S2

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

V

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLAQR2 is identical to CLAQR3 except that it avoids
recursion by calling CLAHQR instead of CLAQR4.
Aggressive early deflation:
This subroutine 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

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

WANTZ

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

N

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

KTOP

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

KBOT

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

NW

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

H

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

LDH

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

ILOZ

ILOZ is INTEGER

IHIZ

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

Z

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

LDZ

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

NS

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

ND

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

SH

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

V

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

LDV

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

NH

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

T

T is COMPLEX array, dimension (LDT,NW)

LDT

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

NV

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

WV

WV is COMPLEX array, dimension (LDWV,NW)

LDWV

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

WORK

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

LWORK

LWORK is INTEGER
The dimension of the work array WORK.  LWORK = 2*NW
suffices, but greater efficiency may result from larger
values of LWORK.
If LWORK = -1, then a workspace query is assumed; CLAQR2
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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

Aggressive early deflation:
CLAQR3 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

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

WANTZ

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

N

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

KTOP

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

KBOT

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

NW

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

H

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

LDH

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

ILOZ

ILOZ is INTEGER

IHIZ

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

Z

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

LDZ

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

NS

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

ND

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

SH

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

V

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

LDV

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

NH

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

T

T is COMPLEX array, dimension (LDT,NW)

LDT

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

NV

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

WV

WV is COMPLEX array, dimension (LDWV,NW)

LDWV

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

WORK

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

LWORK

LWORK is INTEGER
The dimension of the work array WORK.  LWORK = 2*NW
suffices, but greater efficiency may result from larger
values of LWORK.
If LWORK = -1, then a workspace query is assumed; CLAQR3
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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLAQR4 implements one level of recursion for CLAQR0.
It is a complete implementation of the small bulge multi-shift
QR algorithm.  It may be called by CLAQR0 and, for large enough
deflation window size, it may be called by CLAQR3.  This
subroutine is identical to CLAQR0 except that it calls CLAQR2
instead of CLAQR3.
CLAQR4 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

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

WANTZ

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

N

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

ILO

ILO is INTEGER

IHI

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

H

H is COMPLEX array, dimension (LDH,N)
 On entry, the upper Hessenberg matrix H.
 On exit, if INFO = 0 and WANTT is .TRUE., then H
 contains the upper triangular matrix T from the Schur
 decomposition (the Schur form). If INFO = 0 and WANT is
 .FALSE., then the contents of H are unspecified on exit.
 (The output value of H when INFO > 0 is given under the
 description of INFO below.)
 This subroutine may explicitly set H(i,j) = 0 for i > j and
 j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.

LDH

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

W

W is COMPLEX array, dimension (N)
 The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
 in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
 stored in the same order as on the diagonal of the Schur
 form returned in H, with W(i) = H(i,i).

ILOZ

ILOZ is INTEGER

IHIZ

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

Z

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

LDZ

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

WORK

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

LWORK

LWORK is INTEGER
 The dimension of the array WORK.  LWORK >= max(1,N)
 is sufficient, but LWORK typically as large as 6*N may
 be required for optimal performance.  A workspace query
 to determine the optimal workspace size is recommended.
 If LWORK = -1, then CLAQR4 does a workspace query.
 In this case, CLAQR4 checks the input parameters and
 estimates the optimal workspace size for the given
 values of N, ILO and IHI.  The estimate is returned
 in WORK(1).  No error message related to LWORK is
 issued by XERBLA.  Neither H nor Z are accessed.

INFO

INFO is INTEGER
   = 0:  successful exit
   > 0:  if INFO = i, CLAQR4 failed to compute all of
      the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
      and WI contain those eigenvalues which have been
      successfully computed.  (Failures are rare.)
      If INFO > 0 and WANT is .FALSE., then on exit,
      the remaining unconverged eigenvalues are the eigen-
      values of the upper Hessenberg matrix rows and
      columns ILO through INFO of the final, output
      value of H.
      If INFO > 0 and WANTT is .TRUE., then on exit
 (*)  (initial value of H)*U  = U*(final value of H)
      where U is a unitary matrix.  The final
      value of  H is upper Hessenberg and triangular in
      rows and columns INFO+1 through IHI.
      If INFO > 0 and WANTZ is .TRUE., then on exit
        (final value of Z(ILO:IHI,ILOZ:IHIZ)
         =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
      where U is the unitary matrix in (*) (regard-
      less of the value of WANTT.)
      If INFO > 0 and WANTZ is .FALSE., then Z is not
      accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR Algorithm Part II: Aggressive Early Deflation, SIAM Journal of Matrix Analysis, volume 23, pages 948--973, 2002.

CLAQR5 called by CLAQR0 performs a
single small-bulge multi-shift QR sweep.

WANTT

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

WANTZ

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

KACC22

  KACC22 is INTEGER with value 0, 1, or 2.
     Specifies the computation mode of far-from-diagonal
     orthogonal updates.
= 0: CLAQR5 does not accumulate reflections and does not
     use matrix-matrix multiply to update far-from-diagonal
     matrix entries.
= 1: CLAQR5 accumulates reflections and uses matrix-matrix
     multiply to update the far-from-diagonal matrix entries.
= 2: Same as KACC22 = 1. This option used to enable exploiting
     the 2-by-2 structure during matrix multiplications, but
     this is no longer supported.

N

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

KTOP

KTOP is INTEGER

KBOT

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

NSHFTS

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

S

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

H

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

LDH

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

ILOZ

ILOZ is INTEGER

IHIZ

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

Z

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

LDZ

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

V

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

LDV

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

U

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

LDU

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

NV

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

WV

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

LDWV

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

NH

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

WH

WH is COMPLEX array, dimension (LDWH,NH)

LDWH

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

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

UPLO

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

N

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

KD

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

AB

AB is COMPLEX array, dimension (LDAB,N)
On entry, the upper or lower triangle of the symmetric band
matrix A, stored in the first KD+1 rows of the array.  The
j-th column of A is stored in the j-th column of the array AB
as follows:
if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
On exit, if INFO = 0, the triangular factor U or L from the
Cholesky factorization A = U**H *U or A = L*L**H of the band
matrix A, in the same storage format as A.

LDAB

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

S

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

SCOND

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

AMAX

AMAX is REAL
Absolute value of largest matrix entry.

EQUED

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

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

UPLO

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

N

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

AP

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

S

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

SCOND

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

AMAX

AMAX is REAL
Absolute value of largest matrix entry.

EQUED

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

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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAR1V 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

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

B1

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

BN

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

LAMBDA

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

L

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

D

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

LD

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

LLD

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

PIVMIN

PIVMIN is REAL
 The minimum pivot in the Sturm sequence.

GAPTOL

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

Z

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

WANTNC

WANTNC is LOGICAL
 Specifies whether NEGCNT has to be computed.

NEGCNT

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

ZTZ

ZTZ is REAL
 The square of the 2-norm of Z.

MINGMA

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

R

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

ISUPPZ

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

NRMINV

NRMINV is REAL
 NRMINV = 1/SQRT( ZTZ )

RESID

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

RQCORR

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

WORK

WORK is REAL array, dimension (4*N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLAR2V 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

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

X

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

Y

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

Z

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

INCX

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

C

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

S

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

INCC

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARCM 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

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

N

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

A

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

LDA

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

B

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

LDB

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

C

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

LDC

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

RWORK

RWORK is REAL array, dimension (2*M*N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARF 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 (the conjugate transpose of H), supply conjg(tau) instead
tau.

SIDE

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

M

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

N

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

V

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

INCV

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

TAU

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

C

C is COMPLEX array, dimension (LDC,N)
On entry, the M-by-N matrix C.
On exit, C is overwritten by the matrix H * C if SIDE = 'L',
or C * H if SIDE = 'R'.

LDC

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARFB 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

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

TRANS

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

DIRECT

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

STOREV

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

M

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

N

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

K

K is INTEGER
The order of the matrix T (= the number of elementary
reflectors whose product defines the block reflector).
If SIDE = 'L', M >= K >= 0;
if SIDE = 'R', N >= K >= 0.

V

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

LDV

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

T

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

LDT

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

C

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

LDC

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

WORK

WORK is COMPLEX array, dimension (LDWORK,K)

LDWORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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 )

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

IDENT

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

M

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

N

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

K

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

T

T is COMPLEX array, dimension (LDT,K)
The upper-triangular K-by-K matrix T in the representation
of the block reflector.

LDT

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

A

A is COMPLEX array, dimension (LDA,N)
On entry:
 a) In the K-by-N upper-trapezoidal part A: input matrix A.
 b) In the columns below the diagonal: columns of V1
    (ones are not stored on the diagonal).
On exit:
  A is overwritten by rectangular K-by-N product H*A.
See Further Details section.

LDA

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

B

B is COMPLEX array, dimension (LDB,N)
On entry:
  a) In the M-by-(N-K) right block: input matrix B.
  b) In the M-by-N left block: columns of V2.
On exit:
  B is overwritten by rectangular M-by-N product H*B.
See Further Details section.

LDB

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

WORK

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

LDWORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

(1) Description of the Algebraic Operation.
The matrix A is a K-by-N matrix composed of two column block
matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):
A = ( A1, A2 ).
The matrix B is an M-by-N matrix composed of two column block
matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):
B = ( B1, B2 ).
Perform the operation:
   ( A_out ) := H * ( A_in ) = ( I - V * T * V**H ) * ( A_in ) =
   ( B_out )        ( B_in )                          ( B_in )
              = ( I - ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in )
                      ( V2 )                            ( B_in )
 On input:
a) ( A_in )  consists of two block columns:
   ( B_in )
   ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))
   ( B_in )   (( B1_in ) ( B2_in ))   ((     0 ) ( B2_in )),
   where the column blocks are:
   (  A1_in )  is a K-by-K upper-triangular matrix stored in the
               upper triangular part of the array A(1:K,1:K).
   (  B1_in )  is an M-by-K rectangular ZERO matrix and not stored.
   ( A2_in )  is a K-by-(N-K) rectangular matrix stored
              in the array A(1:K,K+1:N).
   ( B2_in )  is an M-by-(N-K) rectangular matrix stored
              in the array B(1:M,K+1:N).
b) V = ( V1 )
       ( V2 )
   where:
   1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;
   2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,
      stored in the lower-triangular part of the array
      A(1:K,1:K) (ones are not stored),
   and V2 is an M-by-K rectangular stored the array B(1:M,1:K),
             (because on input B1_in is a rectangular zero
              matrix that is not stored and the space is
              used to store V2).
c) T is a K-by-K upper-triangular matrix stored
   in the array T(1:K,1:K).
On output:
a) ( A_out ) consists of two  block columns:
   ( B_out )
   ( A_out ) = (( A1_out ) ( A2_out ))
   ( B_out )   (( B1_out ) ( B2_out )),
   where the column blocks are:
   ( A1_out )  is a K-by-K square matrix, or a K-by-K
               upper-triangular matrix, if V1 is an
               identity matrix. AiOut is stored in
               the array A(1:K,1:K).
   ( B1_out )  is an M-by-K rectangular matrix stored
               in the array B(1:M,K:N).
   ( A2_out )  is a K-by-(N-K) rectangular matrix stored
               in the array A(1:K,K+1:N).
   ( B2_out )  is an M-by-(N-K) rectangular matrix stored
               in the array B(1:M,K+1:N).
The operation above can be represented as the same operation
on each block column:
   ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**H ) * ( A1_in )
   ( B1_out )        (     0 )                          (     0 )
   ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**H ) * ( A2_in )
   ( B2_out )        ( B2_in )                          ( B2_in )
If IDENT != 'I':
   The computation for column block 1:
   A1_out: = A1_in - V1*T*(V1**H)*A1_in
   B1_out: = - V2*T*(V1**H)*A1_in
   The computation for column block 2, which exists if N > K:
   A2_out: = A2_in - V1*T*( (V1**H)*A2_in + (V2**H)*B2_in )
   B2_out: = B2_in - V2*T*( (V1**H)*A2_in + (V2**H)*B2_in )
If IDENT == 'I':
   The operation for column block 1:
   A1_out: = A1_in - V1*T*A1_in
   B1_out: = - V2*T*A1_in
   The computation for column block 2, which exists if N > K:
   A2_out: = A2_in - T*( A2_in + (V2**H)*B2_in )
   B2_out: = B2_in - V2*T*( A2_in + (V2**H)*B2_in )
(2) Description of the Algorithmic Computation.
In the first step, we compute column block 2, i.e. A2 and B2.
Here, we need to use the K-by-(N-K) rectangular workspace
matrix W2 that is of the same size as the matrix A2.
W2 is stored in the array WORK(1:K,1:(N-K)).
In the second step, we compute column block 1, i.e. A1 and B1.
Here, we need to use the K-by-K square workspace matrix W1
that is of the same size as the as the matrix A1.
W1 is stored in the array WORK(1:K,1:K).
NOTE: Hence, in this routine, we need the workspace array WORK
only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from
the first step and W1 from the second step.
Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I',
more computations than in the Case (B).
if( IDENT != 'I' ) then
 if ( N > K ) then
   (First Step - column block 2)
   col2_(1) W2: = A2
   col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2
   col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
   col2_(4) W2: = T * W2
   col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
   col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
   col2_(7) A2: = A2 - W2
 else
   (Second Step - column block 1)
   col1_(1) W1: = A1
   col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1
   col1_(3) W1: = T * W1
   col1_(4) B1: = - V2 * W1 = - B1 * W1
   col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
   col1_(6) square A1: = A1 - W1
 end if
end if
Case (B), when V1 is an identity matrix, i.e. IDENT == 'I',
less computations than in the Case (A)
if( IDENT == 'I' ) then
 if ( N > K ) then
   (First Step - column block 2)
   col2_(1) W2: = A2
   col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
   col2_(4) W2: = T * W2
   col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
   col2_(7) A2: = A2 - W2
 else
   (Second Step - column block 1)
   col1_(1) W1: = A1
   col1_(3) W1: = T * W1
   col1_(4) B1: = - V2 * W1 = - B1 * W1
   col1_(6) upper-triangular_of_(A1): = A1 - W1
 end if
end if
Combine these cases (A) and (B) together, this is the resulting
algorithm:
if ( N > K ) then
  (First Step - column block 2)
  col2_(1)  W2: = A2
  if( IDENT != 'I' ) then
    col2_(2)  W2: = (V1**H) * W2
                  = (unit_lower_tr_of_(A1)**H) * W2
  end if
  col2_(3)  W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2]
  col2_(4)  W2: = T * W2
  col2_(5)  B2: = B2 - V2 * W2 = B2 - B1 * W2
  if( IDENT != 'I' ) then
    col2_(6)    W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
  end if
  col2_(7) A2: = A2 - W2
else
(Second Step - column block 1)
  col1_(1) W1: = A1
  if( IDENT != 'I' ) then
    col1_(2) W1: = (V1**H) * W1
                = (unit_lower_tr_of_(A1)**H) * W1
  end if
  col1_(3) W1: = T * W1
  col1_(4) B1: = - V2 * W1 = - B1 * W1
  if( IDENT != 'I' ) then
    col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
    col1_(6_a) below_diag_of_(A1): =  - below_diag_of_(W1)
  end if
  col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)
end if

CLARFG 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

N is INTEGER
The order of the elementary reflector.

ALPHA

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

X

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

INCX

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

TAU

TAU is COMPLEX
The value tau.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARFGP 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

N is INTEGER
The order of the elementary reflector.

ALPHA

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

X

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

INCX

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

TAU

TAU is COMPLEX
The value tau.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARFT 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

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

STOREV

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

N

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

K

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

V

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

LDV

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

TAU

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

T

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

LDT

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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 )

CLARFX 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

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

M

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

N

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

V

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

TAU

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

C

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

LDC

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

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARFY applies an elementary reflector, or Householder matrix, H,
to an n x n Hermitian matrix C, from both the left and the right.
H is represented in the form
   H = I - tau * v * v'
where  tau  is a scalar and  v  is a vector.
If  tau  is  zero, then  H  is taken to be the unit matrix.

UPLO

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

N

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

V

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

INCV

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

TAU

TAU is COMPLEX
The value tau as described above.

C

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

LDC

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

WORK

WORK is COMPLEX array, dimension (N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLARGV 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 CLARTG,
but differ from the BLAS1 routine CROTG):
   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

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

X

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

INCX

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

Y

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

INCY

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

C

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

INCC

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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.

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

IDIST

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

ISEED

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

N

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

X

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

This routine calls the auxiliary routine SLARUV 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.

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

N

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

VL

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

VU

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

D

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

L

L is REAL 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 SLARRE.
On exit, L is overwritten.

PIVMIN

PIVMIN is REAL
The minimum pivot allowed in the Sturm sequence.

ISPLIT

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

M

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

DOL

DOL is INTEGER

DOU

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

MINRGP

MINRGP is REAL

RTOL1

RTOL1 is REAL

RTOL2

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

W

W is REAL 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 SLARRE is expected here ). Furthermore, they are with
respect to the shift of the corresponding root representation
for their block. On exit, W holds the eigenvalues of the
UNshifted matrix.

WERR

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

WGAP

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

IBLOCK

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

INDEXW

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

GERS

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

Z

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

LDZ

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

ISUPPZ

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

WORK

WORK is REAL array, dimension (12*N)

IWORK

IWORK is INTEGER array, dimension (7*N)

INFO

INFO is INTEGER
= 0:  successful exit
> 0:  A problem occurred in CLARRV.
< 0:  One of the called subroutines signaled an internal problem.
      Needs inspection of the corresponding parameter IINFO
      for further information.
=-1:  Problem in SLARRB when refining a child's eigenvalues.
=-2:  Problem in SLARRF 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 SLARRB 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.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLARTV 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

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

X

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

INCX

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

Y

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

INCY

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

C

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

S

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

INCC

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLASCL 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

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 CGBTRF for storage details.

KL

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

KU

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

CFROM

CFROM is REAL

CTO

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

M

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

N

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

A

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

LDA

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

UPLO

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

M

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

N

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

ALPHA

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

BETA

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

A

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

LDA

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLASR 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

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

PIVOT

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

DIRECT

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

M

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

N

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

C

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

S

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

A

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

LDA

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLASWP 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

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

A

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

LDA

LDA is INTEGER
The leading dimension of the array A.

K1

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

K2

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

IPIV

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

INCX

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLATBS 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 CTBSV 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

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

TRANS

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

DIAG

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

NORMIN

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

N

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

KD

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

AB

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

LDAB

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

X

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

SCALE

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

CNORM

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A rough bound on x is computed; if that is less than overflow, CTBSV
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 CTBSV 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 CTBSV if 1/M(n) and 1/G(n) are both greater
than max(underflow, 1/overflow).

CLATDF 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 CGETC2. On entry RHS = f holds the
contribution from earlier solved sub-systems, and on return RHS = x.
The factorization of Z returned by CGETC2 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

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

N

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

Z

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

LDZ

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

RHS

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

RDSUM

RDSUM is REAL
On entry, the sum of squares of computed contributions to
the Dif-estimate under computation by CTGSYL, 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 CTGSY2 is called by CTGSYL.

RDSCAL

RDSCAL is REAL
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 CTGSY2 is called by
CTGSYL.

IPIV

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

JPIV

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

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

CLATPS 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 CTPSV 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

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

TRANS

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

DIAG

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

NORMIN

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

N

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

AP

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

X

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

SCALE

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

CNORM

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A rough bound on x is computed; if that is less than overflow, CTPSV
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 CTPSV 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 CTPSV if 1/M(n) and 1/G(n) are both greater
than max(underflow, 1/overflow).

CLATRD 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', CLATRD reduces the last NB rows and columns of a
matrix, of which the upper triangle is supplied;
if UPLO = 'L', CLATRD reduces the first NB rows and columns of a
matrix, of which the lower triangle is supplied.
This is an auxiliary routine called by CHETRD.

UPLO

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

N

N is INTEGER
The order of the matrix A.

NB

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

A

A is COMPLEX array, dimension (LDA,N)
On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
n-by-n upper triangular part of A contains the upper
triangular part of the matrix A, and the strictly lower
triangular part of A is not referenced.  If UPLO = 'L', the
leading n-by-n lower triangular part of A contains the lower
triangular part of the matrix A, and the strictly upper
triangular part of A is not referenced.
On exit:
if UPLO = 'U', the last NB columns have been reduced to
  tridiagonal form, with the diagonal elements overwriting
  the diagonal elements of A; the elements above the diagonal
  with the array TAU, represent the unitary matrix Q as a
  product of elementary reflectors;
if UPLO = 'L', the first NB columns have been reduced to
  tridiagonal form, with the diagonal elements overwriting
  the diagonal elements of A; the elements below the diagonal
  with the array TAU, represent the  unitary matrix Q as a
  product of elementary reflectors.
See Further Details.

LDA

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

E

E is REAL array, dimension (N-1)
If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
elements of the last NB columns of the reduced matrix;
if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
the first NB columns of the reduced matrix.

TAU

TAU is COMPLEX array, dimension (N-1)
The scalar factors of the elementary reflectors, stored in
TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
See Further Details.

W

W is COMPLEX array, dimension (LDW,NB)
The n-by-nb matrix W required to update the unreduced part
of A.

LDW

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

CLATRS 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
CTRSV 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

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

TRANS

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

DIAG

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

NORMIN

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

N

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

A

A is COMPLEX array, dimension (LDA,N)
The triangular matrix A.  If UPLO = 'U', the leading n by n
upper triangular part of the array A contains the upper
triangular matrix, and the strictly lower triangular part of
A is not referenced.  If UPLO = 'L', the leading n by n lower
triangular part of the array A contains the lower triangular
matrix, and the strictly upper triangular part of A is not
referenced.  If DIAG = 'U', the diagonal elements of A are
also not referenced and are assumed to be 1.

LDA

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

X

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

SCALE

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

CNORM

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A rough bound on x is computed; if that is less than overflow, CTRSV
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 CTRSV 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 CTRSV if 1/M(n) and 1/G(n) are both greater
than max(underflow, 1/overflow).

CLAUU2 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

UPLO is CHARACTER*1
Specifies whether the triangular factor stored in the array A
is upper or lower triangular:
= 'U':  Upper triangular
= 'L':  Lower triangular

N

N is INTEGER
The order of the triangular factor U or L.  N >= 0.

A

A is COMPLEX array, dimension (LDA,N)
On entry, the triangular factor U or L.
On exit, if UPLO = 'U', the upper triangle of A is
overwritten with the upper triangle of the product U * U**H;
if UPLO = 'L', the lower triangle of A is overwritten with
the lower triangle of the product L**H * L.

LDA

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CLAUUM 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

UPLO is CHARACTER*1
Specifies whether the triangular factor stored in the array A
is upper or lower triangular:
= 'U':  Upper triangular
= 'L':  Lower triangular

N

N is INTEGER
The order of the triangular factor U or L.  N >= 0.

A

A is COMPLEX array, dimension (LDA,N)
On entry, the triangular factor U or L.
On exit, if UPLO = 'U', the upper triangle of A is
overwritten with the upper triangle of the product U * U**H;
if UPLO = 'L', the lower triangle of A is overwritten with
the lower triangle of the product L**H * L.

LDA

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

INFO

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CROT   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

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

CX

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

INCX

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

CY

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

INCY

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

C

C is REAL

S

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CSPMV  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

UPLO is CHARACTER*1
 On entry, UPLO specifies whether the upper or lower
 triangular part of the matrix A is supplied in the packed
 array AP as follows:
    UPLO = 'U' or 'u'   The upper triangular part of A is
                        supplied in AP.
    UPLO = 'L' or 'l'   The lower triangular part of A is
                        supplied in AP.
 Unchanged on exit.

N

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

ALPHA

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

AP

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

X

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

INCX

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

BETA

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

Y

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

INCY

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CSPR    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

UPLO is CHARACTER*1
 On entry, UPLO specifies whether the upper or lower
 triangular part of the matrix A is supplied in the packed
 array AP as follows:
    UPLO = 'U' or 'u'   The upper triangular part of A is
                        supplied in AP.
    UPLO = 'L' or 'l'   The lower triangular part of A is
                        supplied in AP.
 Unchanged on exit.

N

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

ALPHA

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

X

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

INCX

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

AP

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CSRSCL 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

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

SA

SA is REAL
The scalar a which is used to divide each component of x.
SA must be >= 0, or the subroutine will divide by zero.

SX

SX is COMPLEX array, dimension
               (1+(N-1)*abs(INCX))
The n-element vector x.

INCX

INCX is INTEGER
The increment between successive values of the vector SX.
> 0:  SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i),     1< i<= n

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

CTPRFB 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

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

TRANS

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

DIRECT

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

STOREV

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

M

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

N

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

K

K is INTEGER
The order of the matrix T, i.e. the number of elementary
reflectors whose product defines the block reflector.
K >= 0.

L

L is INTEGER
The order of the trapezoidal part of V.
K >= L >= 0.  See Further Details.

V

V is COMPLEX array, dimension
                      (LDV,K) if STOREV = 'C'
                      (LDV,M) if STOREV = 'R' and SIDE = 'L'
                      (LDV,N) if STOREV = 'R' and SIDE = 'R'
The pentagonal matrix V, which contains the elementary reflectors
H(1), H(2), ..., H(K).  See Further Details.

LDV

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

T

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

LDT

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

A

A is COMPLEX array, dimension
(LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
On entry, the K-by-N or M-by-K matrix A.
On exit, A is overwritten by the corresponding block of
H*C or H**H*C or C*H or C*H**H.  See Further Details.

LDA

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

B

B is COMPLEX array, dimension (LDB,N)
On entry, the M-by-N matrix B.
On exit, B is overwritten by the corresponding block of
H*C or H**H*C or C*H or C*H**H.  See Further Details.

LDB

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

WORK

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

LDWORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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.

ICMAX1 finds the index of the first vector element of maximum absolute value.
Based on ICAMAX from Level 1 BLAS.
The change is to use the 'genuine' absolute value.

N

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

CX

CX is COMPLEX array, dimension (N)
The vector CX. The ICMAX1 function returns the index of its first
element of maximum absolute value.

INCX

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Nick Higham for use with CLACON.

ILACLC scans A for its last non-zero column.

M

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

N

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

A

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

LDA

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

ILACLR scans A for its last non-zero row.

M

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

N

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

A

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

LDA

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

IZMAX1 finds the index of the first vector element of maximum absolute value.
Based on IZAMAX from Level 1 BLAS.
The change is to use the 'genuine' absolute value.

N

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

ZX

ZX is COMPLEX*16 array, dimension (N)
The vector ZX. The IZMAX1 function returns the index of its first
element of maximum absolute value.

INCX

INCX is INTEGER
The spacing between successive values of ZX.  INCX >= 1.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Nick Higham for use with ZLACON.

SCSUM1 takes the sum of the absolute values of a complex
vector and returns a single precision result.
Based on SCASUM from the Level 1 BLAS.
The change is to use the 'genuine' absolute value.

N

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

CX

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

INCX

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Nick Higham for use with CLACON.