## table of contents

complexOTHERauxiliary(3) | LAPACK | complexOTHERauxiliary(3) |

# NAME¶

complexOTHERauxiliary - complex

# SYNOPSIS¶

## Functions¶

subroutine **clabrd** (M, N, NB, A, LDA, D, E, TAUQ, TAUP, X,
LDX, Y, LDY)

**CLABRD** reduces the first nb rows and columns of a general matrix to a
bidiagonal form. subroutine **clacgv** (N, X, INCX)

**CLACGV** conjugates a complex vector. subroutine **clacn2** (N, V, X,
EST, KASE, ISAVE)

**CLACN2** estimates the 1-norm of a square matrix, using reverse
communication for evaluating matrix-vector products. subroutine
**clacon** (N, V, X, EST, KASE)

**CLACON** estimates the 1-norm of a square matrix, using reverse
communication for evaluating matrix-vector products. subroutine
**clacp2** (UPLO, M, N, A, LDA, B, LDB)

**CLACP2** copies all or part of a real two-dimensional array to a complex
array. subroutine **clacpy** (UPLO, M, N, A, LDA, B, LDB)

**CLACPY** copies all or part of one two-dimensional array to another.
subroutine **clacrm** (M, N, A, LDA, B, LDB, C, LDC, RWORK)

**CLACRM** multiplies a complex matrix by a square real matrix. subroutine
**clacrt** (N, CX, INCX, CY, INCY, C, S)

**CLACRT** performs a linear transformation of a pair of complex vectors.
complex function **cladiv** (X, Y)

**CLADIV** performs complex division in real arithmetic, avoiding
unnecessary overflow. subroutine **claein** (RIGHTV, NOINIT, N, H, LDH,
W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO)

**CLAEIN** computes a specified right or left eigenvector of an upper
Hessenberg matrix by inverse iteration. subroutine **claev2** (A, B, C,
RT1, RT2, CS1, SN1)

**CLAEV2** computes the eigenvalues and eigenvectors of a 2-by-2
symmetric/Hermitian matrix. subroutine **clags2** (UPPER, A1, A2, A3, B1,
B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)

**CLAGS2** subroutine **clagtm** (TRANS, N, NRHS, ALPHA, DL, D, DU, X,
LDX, BETA, B, LDB)

**CLAGTM** performs a matrix-matrix product of the form C =
αAB+βC, where A is a tridiagonal matrix, B and C are
rectangular matrices, and α and β are scalars, which may be 0,
1, or -1. subroutine **clahqr** (WANTT, WANTZ, N, ILO, IHI, H, LDH, W,
ILOZ, IHIZ, Z, LDZ, INFO)

**CLAHQR** computes the eigenvalues and Schur factorization of an upper
Hessenberg matrix, using the double-shift/single-shift QR algorithm.
subroutine **clahr2** (N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)

**CLAHR2** reduces the specified number of first columns of a general
rectangular matrix A so that elements below the specified subdiagonal are
zero, and returns auxiliary matrices which are needed to apply the
transformation to the unreduced part of A. subroutine **claic1** (JOB, J,
X, SEST, W, GAMMA, SESTPR, S, C)

**CLAIC1** applies one step of incremental condition estimation. real
function **clangt** (NORM, N, DL, D, DU)

**CLANGT** returns the value of the 1-norm, Frobenius norm, infinity-norm,
or the largest absolute value of any element of a general tridiagonal
matrix. real function **clanhb** (NORM, UPLO, N, K, AB, LDAB, WORK)

**CLANHB** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a Hermitian band
matrix. real function **clanhp** (NORM, UPLO, N, AP, WORK)

**CLANHP** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a complex
Hermitian matrix supplied in packed form. real function **clanhs** (NORM,
N, A, LDA, WORK)

**CLANHS** returns the value of the 1-norm, Frobenius norm, infinity-norm,
or the largest absolute value of any element of an upper Hessenberg matrix.
real function **clanht** (NORM, N, D, E)

**CLANHT** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a complex
Hermitian tridiagonal matrix. real function **clansb** (NORM, UPLO, N, K,
AB, LDAB, WORK)

**CLANSB** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a symmetric band
matrix. real function **clansp** (NORM, UPLO, N, AP, WORK)

**CLANSP** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a symmetric
matrix supplied in packed form. real function **clantb** (NORM, UPLO,
DIAG, N, K, AB, LDAB, WORK)

**CLANTB** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a triangular band
matrix. real function **clantp** (NORM, UPLO, DIAG, N, AP, WORK)

**CLANTP** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a triangular
matrix supplied in packed form. real function **clantr** (NORM, UPLO,
DIAG, M, N, A, LDA, WORK)

**CLANTR** returns the value of the 1-norm, or the Frobenius norm, or the
infinity norm, or the element of largest absolute value of a trapezoidal or
triangular matrix. subroutine **clapll** (N, X, INCX, Y, INCY, SSMIN)

**CLAPLL** measures the linear dependence of two vectors. subroutine
**clapmr** (FORWRD, M, N, X, LDX, K)

**CLAPMR** rearranges rows of a matrix as specified by a permutation
vector. subroutine **clapmt** (FORWRD, M, N, X, LDX, K)

**CLAPMT** performs a forward or backward permutation of the columns of a
matrix. subroutine **claqhb** (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX,
EQUED)

**CLAQHB** scales a Hermitian band matrix, using scaling factors computed
by cpbequ. subroutine **claqhp** (UPLO, N, AP, S, SCOND, AMAX, EQUED)

**CLAQHP** scales a Hermitian matrix stored in packed form. subroutine
**claqp2** (M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)

**CLAQP2** computes a QR factorization with column pivoting of the matrix
block. subroutine **claqps** (M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU,
VN1, VN2, AUXV, F, LDF)

**CLAQPS** computes a step of QR factorization with column pivoting of a
real m-by-n matrix A by using BLAS level 3. subroutine **claqr0** (WANTT,
WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)

**CLAQR0** computes the eigenvalues of a Hessenberg matrix, and optionally
the matrices from the Schur decomposition. subroutine **claqr1** (N, H,
LDH, S1, S2, V)

**CLAQR1** sets a scalar multiple of the first column of the product of
2-by-2 or 3-by-3 matrix H and specified shifts. subroutine **claqr2**
(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V,
LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)

**CLAQR2** performs the unitary similarity transformation of a Hessenberg
matrix to detect and deflate fully converged eigenvalues from a trailing
principal submatrix (aggressive early deflation). subroutine **claqr3**
(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V,
LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)

**CLAQR3** performs the unitary similarity transformation of a Hessenberg
matrix to detect and deflate fully converged eigenvalues from a trailing
principal submatrix (aggressive early deflation). subroutine **claqr4**
(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK,
INFO)

**CLAQR4** computes the eigenvalues of a Hessenberg matrix, and optionally
the matrices from the Schur decomposition. subroutine **claqr5** (WANTT,
WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV,
U, LDU, NV, WV, LDWV, NH, WH, LDWH)

**CLAQR5** performs a single small-bulge multi-shift QR sweep. subroutine
**claqsb** (UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)

**CLAQSB** scales a symmetric/Hermitian band matrix, using scaling factors
computed by spbequ. subroutine **claqsp** (UPLO, N, AP, S, SCOND, AMAX,
EQUED)

**CLAQSP** scales a symmetric/Hermitian matrix in packed storage, using
scaling factors computed by sppequ. subroutine **clar1v** (N, B1, BN,
LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R,
ISUPPZ, NRMINV, RESID, RQCORR, WORK)

**CLAR1V** computes the (scaled) r-th column of the inverse of the
submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.
subroutine **clar2v** (N, X, Y, Z, INCX, C, S, INCC)

**CLAR2V** applies a vector of plane rotations with real cosines and
complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian
matrices. subroutine **clarcm** (M, N, A, LDA, B, LDB, C, LDC, RWORK)

**CLARCM** copies all or part of a real two-dimensional array to a complex
array. subroutine **clarf** (SIDE, M, N, V, INCV, TAU, C, LDC, WORK)

**CLARF** applies an elementary reflector to a general rectangular matrix.
subroutine **clarfb** (SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T,
LDT, C, LDC, WORK, LDWORK)

**CLARFB** applies a block reflector or its conjugate-transpose to a
general rectangular matrix. subroutine **clarfb_gett** (IDENT, M, N, K,
T, LDT, A, LDA, B, LDB, WORK, LDWORK)

**CLARFB_GETT** subroutine **clarfg** (N, ALPHA, X, INCX, TAU)

**CLARFG** generates an elementary reflector (Householder matrix).
subroutine **clarfgp** (N, ALPHA, X, INCX, TAU)

**CLARFGP** generates an elementary reflector (Householder matrix) with
non-negative beta. subroutine **clarft** (DIRECT, STOREV, N, K, V, LDV,
TAU, T, LDT)

**CLARFT** forms the triangular factor T of a block reflector H = I - vtvH
subroutine **clarfx** (SIDE, M, N, V, TAU, C, LDC, WORK)

**CLARFX** applies an elementary reflector to a general rectangular matrix,
with loop unrolling when the reflector has order ≤ 10. subroutine
**clarfy** (UPLO, N, V, INCV, TAU, C, LDC, WORK)

**CLARFY** subroutine **clargv** (N, X, INCX, Y, INCY, C, INCC)

**CLARGV** generates a vector of plane rotations with real cosines and
complex sines. subroutine **clarnv** (IDIST, ISEED, N, X)

**CLARNV** returns a vector of random numbers from a uniform or normal
distribution. subroutine **clarrv** (N, VL, VU, D, L, PIVMIN, ISPLIT, M,
DOL, DOU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ,
ISUPPZ, WORK, IWORK, INFO)

**CLARRV** computes the eigenvectors of the tridiagonal matrix T = L D LT
given L, D and the eigenvalues of L D LT. subroutine **clartv** (N, X,
INCX, Y, INCY, C, S, INCC)

**CLARTV** applies a vector of plane rotations with real cosines and
complex sines to the elements of a pair of vectors. subroutine **clascl**
(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)

**CLASCL** multiplies a general rectangular matrix by a real scalar defined
as cto/cfrom. subroutine **claset** (UPLO, M, N, ALPHA, BETA, A, LDA)

**CLASET** initializes the off-diagonal elements and the diagonal elements
of a matrix to given values. subroutine **clasr** (SIDE, PIVOT, DIRECT,
M, N, C, S, A, LDA)

**CLASR** applies a sequence of plane rotations to a general rectangular
matrix. subroutine **claswp** (N, A, LDA, K1, K2, IPIV, INCX)

**CLASWP** performs a series of row interchanges on a general rectangular
matrix. subroutine **clatbs** (UPLO, TRANS, DIAG, NORMIN, N, KD, AB,
LDAB, X, SCALE, CNORM, INFO)

**CLATBS** solves a triangular banded system of equations. subroutine
**clatdf** (IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)

**CLATDF** uses the LU factorization of the n-by-n matrix computed by
sgetc2 and computes a contribution to the reciprocal Dif-estimate.
subroutine **clatps** (UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM,
INFO)

**CLATPS** solves a triangular system of equations with the matrix held in
packed storage. subroutine **clatrd** (UPLO, N, NB, A, LDA, E, TAU, W,
LDW)

**CLATRD** reduces the first nb rows and columns of a symmetric/Hermitian
matrix A to real tridiagonal form by an unitary similarity transformation.
subroutine **clatrs** (UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
CNORM, INFO)

**CLATRS** solves a triangular system of equations with the scale factor
set to prevent overflow. subroutine **clauu2** (UPLO, N, A, LDA, INFO)

**CLAUU2** computes the product UUH or LHL, where U and L are upper or
lower triangular matrices (unblocked algorithm). subroutine **clauum**
(UPLO, N, A, LDA, INFO)

**CLAUUM** computes the product UUH or LHL, where U and L are upper or
lower triangular matrices (blocked algorithm). subroutine **crot** (N,
CX, INCX, CY, INCY, C, S)

**CROT** applies a plane rotation with real cosine and complex sine to a
pair of complex vectors. subroutine **cspmv** (UPLO, N, ALPHA, AP, X,
INCX, BETA, Y, INCY)

**CSPMV** computes a matrix-vector product for complex vectors using a
complex symmetric packed matrix subroutine **cspr** (UPLO, N, ALPHA, X,
INCX, AP)

**CSPR** performs the symmetrical rank-1 update of a complex symmetric
packed matrix. subroutine **csrscl** (N, SA, SX, INCX)

**CSRSCL** multiplies a vector by the reciprocal of a real scalar.
subroutine **ctprfb** (SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV,
T, LDT, A, LDA, B, LDB, WORK, LDWORK)

**CTPRFB** applies a real or complex 'triangular-pentagonal' blocked
reflector to a real or complex matrix, which is composed of two blocks.
integer function **icmax1** (N, CX, INCX)

**ICMAX1** finds the index of the first vector element of maximum absolute
value. integer function **ilaclc** (M, N, A, LDA)

**ILACLC** scans a matrix for its last non-zero column. integer function
**ilaclr** (M, N, A, LDA)

**ILACLR** scans a matrix for its last non-zero row. integer function
**izmax1** (N, ZX, INCX)

**IZMAX1** finds the index of the first vector element of maximum absolute
value. real function **scsum1** (N, CX, INCX)

**SCSUM1** forms the 1-norm of the complex vector using the true absolute
value.

# Detailed Description¶

This is the group of complex other auxiliary routines

# Function Documentation¶

## subroutine clabrd (integer M, integer N, integer NB, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, real, dimension( * ) E, complex, dimension( * ) TAUQ, complex, dimension( * ) TAUP, complex, dimension( ldx, * ) X, integer LDX, complex, dimension( ldy, * ) Y, integer LDY)¶

**CLABRD** reduces the first nb rows and columns of a general
matrix to a bidiagonal form.

**Purpose:**

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

**Parameters**

*M*

M is INTEGER

The number of rows in the matrix A.

*N*

N is INTEGER

The number of columns in the matrix A.

*NB*

NB is INTEGER

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

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

The matrices Q and P are represented as products of elementary

reflectors:

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

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

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

where tauq and taup are complex scalars, and v and u are complex

vectors.

If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in

A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in

A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in

A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in

A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).

The elements of the vectors v and u together form the m-by-nb matrix

V and the nb-by-n matrix U**H which are needed, with X and Y, to apply

the transformation to the unreduced part of the matrix, using a block

update of the form: A := A - V*Y**H - X*U**H.

The contents of A on exit are illustrated by the following examples

with nb = 2:

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

( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 )

( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 )

( v1 v2 a a a ) ( v1 1 a a a a )

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

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

( v1 v2 a a a )

where a denotes an element of the original matrix which is unchanged,

vi denotes an element of the vector defining H(i), and ui an element

of the vector defining G(i).

## subroutine clacgv (integer N, complex, dimension( * ) X, integer INCX)¶

**CLACGV** conjugates a complex vector.

**Purpose:**

CLACGV conjugates a complex vector of length N.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clacn2 (integer N, complex, dimension( * ) V, complex, dimension( * ) X, real EST, integer KASE, integer, dimension( 3 ) ISAVE)¶

**CLACN2** estimates the 1-norm of a square matrix, using
reverse communication for evaluating matrix-vector products.

**Purpose:**

CLACN2 estimates the 1-norm of a square, complex matrix A.

Reverse communication is used for evaluating matrix-vector products.

**Parameters**

*N*

N is INTEGER

The order of the matrix. N >= 1.

*V*

V is COMPLEX 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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

Originally named CONEST, dated March 16, 1988.

Last modified: April, 1999

This is a thread safe version of 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)

**Contributors:**

**References:**

a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.

## subroutine clacon (integer N, complex, dimension( n ) V, complex, dimension( n ) X, real EST, integer KASE)¶

**CLACON** estimates the 1-norm of a square matrix, using
reverse communication for evaluating matrix-vector products.

**Purpose:**

CLACON estimates the 1-norm of a square, complex matrix A.

Reverse communication is used for evaluating matrix-vector products.

**Parameters**

*N*

N is INTEGER

The order of the matrix. N >= 1.

*V*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

Last modified: April, 1999

**Contributors:**

**References:**

a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.

## subroutine clacp2 (character UPLO, integer M, integer N, real, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)¶

**CLACP2** copies all or part of a real two-dimensional array
to a complex array.

**Purpose:**

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

complex matrix B.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies the part of the matrix A to be copied to B.

= 'U': Upper triangular part

= 'L': Lower triangular part

Otherwise: All of the matrix A

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clacpy (character UPLO, integer M, integer N, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB)¶

**CLACPY** copies all or part of one two-dimensional array to
another.

**Purpose:**

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

matrix B.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies the part of the matrix A to be copied to B.

= 'U': Upper triangular part

= 'L': Lower triangular part

Otherwise: All of the matrix A

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clacrm (integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, complex, dimension( ldc, * ) C, integer LDC, real, dimension( * ) RWORK)¶

**CLACRM** multiplies a complex matrix by a square real
matrix.

**Purpose:**

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.

**Parameters**

*M*

M is INTEGER

The number of rows of the matrix A and of the matrix C.

M >= 0.

*N*

N is INTEGER

The number of columns and rows of the matrix B and

the number of columns of the matrix C.

N >= 0.

*A*

A is COMPLEX 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)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clacrt (integer N, complex, dimension( * ) CX, integer INCX, complex, dimension( * ) CY, integer INCY, complex C, complex S)¶

**CLACRT** performs a linear transformation of a pair of
complex vectors.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## complex function cladiv (complex X, complex Y)¶

**CLADIV** performs complex division in real arithmetic,
avoiding unnecessary overflow.

**Purpose:**

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.

**Parameters**

*X*

X is COMPLEX

*Y*

Y is COMPLEX

The complex scalars X and Y.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claein (logical RIGHTV, logical NOINIT, integer N, complex, dimension( ldh, * ) H, integer LDH, complex W, complex, dimension( * ) V, complex, dimension( ldb, * ) B, integer LDB, real, dimension( * ) RWORK, real EPS3, real SMLNUM, integer INFO)¶

**CLAEIN** computes a specified right or left eigenvector of an
upper Hessenberg matrix by inverse iteration.

**Purpose:**

CLAEIN uses inverse iteration to find a right or left eigenvector

corresponding to the eigenvalue W of a complex upper Hessenberg

matrix H.

**Parameters**

*RIGHTV*

RIGHTV is LOGICAL

= .TRUE. : compute right eigenvector;

= .FALSE.: compute left eigenvector.

*NOINIT*

NOINIT is LOGICAL

= .TRUE. : no initial vector supplied in V

= .FALSE.: initial vector supplied in V.

*N*

N is INTEGER

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

*H*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claev2 (complex A, complex B, complex C, real RT1, real RT2, real CS1, complex SN1)¶

**CLAEV2** computes the eigenvalues and eigenvectors of a
2-by-2 symmetric/Hermitian matrix.

**Purpose:**

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

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

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

RT2 may be inaccurate if there is massive cancellation in the

determinant A*C-B*B; higher precision or correctly rounded or

correctly truncated arithmetic would be needed to compute RT2

accurately in all cases.

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

Overflow is possible only if RT1 is within a factor of 5 of overflow.

Underflow is harmless if the input data is 0 or exceeds

underflow_threshold / macheps.

## subroutine clags2 (logical UPPER, real A1, complex A2, real A3, real B1, complex B2, real B3, real CSU, complex SNU, real CSV, complex SNV, real CSQ, complex SNQ)¶

**CLAGS2**

**Purpose:**

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.

**Parameters**

*UPPER*

UPPER is LOGICAL

= .TRUE.: the input matrices A and B are upper triangular.

= .FALSE.: the input matrices A and B are lower triangular.

*A1*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clagtm (character TRANS, integer N, integer NRHS, real ALPHA, complex, dimension( * ) DL, complex, dimension( * ) D, complex, dimension( * ) DU, complex, dimension( ldx, * ) X, integer LDX, real BETA, complex, dimension( ldb, * ) B, integer LDB)¶

**CLAGTM** performs a matrix-matrix product of the form C =
αAB+βC, where A is a tridiagonal matrix, B and C are
rectangular matrices, and α and β are scalars, which may be 0,
1, or -1.

**Purpose:**

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.

**Parameters**

*TRANS*

TRANS is CHARACTER*1

Specifies the operation applied to A.

= 'N': No transpose, B := alpha * A * X + beta * B

= 'T': Transpose, B := alpha * A**T * X + beta * B

= 'C': Conjugate transpose, B := alpha * A**H * X + beta * B

*N*

N is INTEGER

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

*NRHS*

NRHS is INTEGER

The number of right hand sides, i.e., the number of columns

of the matrices X and B.

*ALPHA*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clahqr (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, complex, dimension( ldh, * ) H, integer LDH, complex, dimension( * ) W, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, integer INFO)¶

**CLAHQR** computes the eigenvalues and Schur factorization of
an upper Hessenberg matrix, using the double-shift/single-shift QR
algorithm.

**Purpose:**

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.

**Parameters**

*WANTT*

WANTT is LOGICAL

= .TRUE. : the full Schur form T is required;

= .FALSE.: only eigenvalues are required.

*WANTZ*

WANTZ is LOGICAL

= .TRUE. : the matrix of Schur vectors Z is required;

= .FALSE.: Schur vectors are not required.

*N*

N is INTEGER

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

*ILO*

ILO is INTEGER

*IHI*

IHI is INTEGER

It is assumed that H is already upper triangular in rows and

columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

02-96 Based on modifications by

David Day, Sandia National Laboratory, USA

12-04 Further modifications by

Ralph Byers, University of Kansas, USA

This is a modified version of 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).

## subroutine clahr2 (integer N, integer K, integer NB, complex, dimension( lda, * ) A, integer LDA, complex, dimension( nb ) TAU, complex, dimension( ldt, nb ) T, integer LDT, complex, dimension( ldy, nb ) Y, integer LDY)¶

**CLAHR2** reduces the specified number of first columns of a
general rectangular matrix A so that elements below the specified
subdiagonal are zero, and returns auxiliary matrices which are needed to
apply the transformation to the unreduced part of A.

**Purpose:**

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.

**Parameters**

*N*

N is INTEGER

The order of the matrix A.

*K*

K is INTEGER

The offset for the reduction. Elements below the k-th

subdiagonal in the first NB columns are reduced to zero.

K < N.

*NB*

NB is INTEGER

The number of columns to be reduced.

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

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

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

Each H(i) has the form

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

where tau is a complex scalar, and v is a complex vector with

v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in

A(i+k+1:n,i), and tau in TAU(i).

The elements of the vectors v together form the (n-k+1)-by-nb matrix

V which is needed, with T and Y, to apply the transformation to the

unreduced part of the matrix, using an update of the form:

A := (I - V*T*V**H) * (A - Y*V**H).

The contents of A on exit are illustrated by the following example

with n = 7, k = 3 and nb = 2:

( a a a a a )

( a a a a a )

( a a a a a )

( h h a a a )

( v1 h a a a )

( v1 v2 a a a )

( v1 v2 a a a )

where a denotes an element of the original matrix A, h denotes a

modified element of the upper Hessenberg matrix H, and vi denotes an

element of the vector defining H(i).

This subroutine is a slight modification of LAPACK-3.0's 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.)

**References:**

performance of reduction to Hessenberg form,' ACM Transactions on Mathematical Software, 32(2):180-194, June 2006.

## subroutine claic1 (integer JOB, integer J, complex, dimension( j ) X, real SEST, complex, dimension( j ) W, complex GAMMA, real SESTPR, complex S, complex C)¶

**CLAIC1** applies one step of incremental condition
estimation.

**Purpose:**

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.

**Parameters**

*JOB*

JOB is INTEGER

= 1: an estimate for the largest singular value is computed.

= 2: an estimate for the smallest singular value is computed.

*J*

J is INTEGER

Length of X and W

*X*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clangt (character NORM, integer N, complex, dimension( * ) DL, complex, dimension( * ) D, complex, dimension( * ) DU)¶

**CLANGT** returns the value of the 1-norm, Frobenius norm,
infinity-norm, or the largest absolute value of any element of a general
tridiagonal matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clanhb (character NORM, character UPLO, integer N, integer K, complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) WORK)¶

**CLANHB** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
Hermitian band matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clanhp (character NORM, character UPLO, integer N, complex, dimension( * ) AP, real, dimension( * ) WORK)¶

**CLANHP** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
complex Hermitian matrix supplied in packed form.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clanhs (character NORM, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) WORK)¶

**CLANHS** returns the value of the 1-norm, Frobenius norm,
infinity-norm, or the largest absolute value of any element of an upper
Hessenberg matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clanht (character NORM, integer N, real, dimension( * ) D, complex, dimension( * ) E)¶

**CLANHT** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
complex Hermitian tridiagonal matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clansb (character NORM, character UPLO, integer N, integer K, complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) WORK)¶

**CLANSB** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
symmetric band matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clansp (character NORM, character UPLO, integer N, complex, dimension( * ) AP, real, dimension( * ) WORK)¶

**CLANSP** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
symmetric matrix supplied in packed form.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clantb (character NORM, character UPLO, character DIAG, integer N, integer K, complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) WORK)¶

**CLANTB** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
triangular band matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clantp (character NORM, character UPLO, character DIAG, integer N, complex, dimension( * ) AP, real, dimension( * ) WORK)¶

**CLANTP** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
triangular matrix supplied in packed form.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## real function clantr (character NORM, character UPLO, character DIAG, integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) WORK)¶

**CLANTR** returns the value of the 1-norm, or the Frobenius
norm, or the infinity norm, or the element of largest absolute value of a
trapezoidal or triangular matrix.

**Purpose:**

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.

**Returns**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clapll (integer N, complex, dimension( * ) X, integer INCX, complex, dimension( * ) Y, integer INCY, real SSMIN)¶

**CLAPLL** measures the linear dependence of two vectors.

**Purpose:**

Given two column vectors X and Y, let

A = ( X Y ).

The subroutine first computes the QR factorization of A = Q*R,

and then computes the SVD of the 2-by-2 upper triangular matrix R.

The smaller singular value of R is returned in SSMIN, which is used

as the measurement of the linear dependency of the vectors X and Y.

**Parameters**

*N*

N is INTEGER

The length of the vectors X and Y.

*X*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clapmr (logical FORWRD, integer M, integer N, complex, dimension( ldx, * ) X, integer LDX, integer, dimension( * ) K)¶

**CLAPMR** rearranges rows of a matrix as specified by a
permutation vector.

**Purpose:**

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.

**Parameters**

*FORWRD*

FORWRD is LOGICAL

= .TRUE., forward permutation

= .FALSE., backward permutation

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*X*

X is COMPLEX array, dimension (LDX,N)

On entry, the M by N matrix X.

On exit, X contains the permuted matrix X.

*LDX*

LDX is INTEGER

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

*K*

K is INTEGER array, dimension (M)

On entry, K contains the permutation vector. K is used as

internal workspace, but reset to its original value on

output.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clapmt (logical FORWRD, integer M, integer N, complex, dimension( ldx, * ) X, integer LDX, integer, dimension( * ) K)¶

**CLAPMT** performs a forward or backward permutation of the
columns of a matrix.

**Purpose:**

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.

**Parameters**

*FORWRD*

FORWRD is LOGICAL

= .TRUE., forward permutation

= .FALSE., backward permutation

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*X*

X is COMPLEX array, dimension (LDX,N)

On entry, the M by N matrix X.

On exit, X contains the permuted matrix X.

*LDX*

LDX is INTEGER

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

*K*

K is INTEGER array, dimension (N)

On entry, K contains the permutation vector. K is used as

internal workspace, but reset to its original value on

output.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claqhb (character UPLO, integer N, integer KD, complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) S, real SCOND, real AMAX, character EQUED)¶

**CLAQHB** scales a Hermitian band matrix, using scaling
factors computed by cpbequ.

**Purpose:**

CLAQHB equilibrates an Hermitian band matrix A using the scaling

factors in the vector S.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the

symmetric matrix A is stored.

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

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

*KD*

KD is INTEGER

The number of super-diagonals of the matrix A if UPLO = 'U',

or the number of sub-diagonals if UPLO = 'L'. KD >= 0.

*AB*

AB is COMPLEX 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).

**Internal Parameters:**

THRESH is a threshold value used to decide if scaling should be done

based on the ratio of the scaling factors. If SCOND < THRESH,

scaling is done.

LARGE and SMALL are threshold values used to decide if scaling should

be done based on the absolute size of the largest matrix element.

If AMAX > LARGE or AMAX < SMALL, scaling is done.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claqhp (character UPLO, integer N, complex, dimension( * ) AP, real, dimension( * ) S, real SCOND, real AMAX, character EQUED)¶

**CLAQHP** scales a Hermitian matrix stored in packed form.

**Purpose:**

CLAQHP equilibrates a Hermitian matrix A using the scaling factors

in the vector S.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the

Hermitian matrix A is stored.

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

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

*AP*

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

**Internal Parameters:**

THRESH is a threshold value used to decide if scaling should be done

based on the ratio of the scaling factors. If SCOND < THRESH,

scaling is done.

LARGE and SMALL are threshold values used to decide if scaling should

be done based on the absolute size of the largest matrix element.

If AMAX > LARGE or AMAX < SMALL, scaling is done.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claqp2 (integer M, integer N, integer OFFSET, complex, dimension( lda, * ) A, integer LDA, integer, dimension( * ) JPVT, complex, dimension( * ) TAU, real, dimension( * ) VN1, real, dimension( * ) VN2, complex, dimension( * ) WORK)¶

**CLAQP2** computes a QR factorization with column pivoting of
the matrix block.

**Purpose:**

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.

**Parameters**

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*OFFSET*

OFFSET is INTEGER

The number of rows of the matrix A that must be pivoted

but no factorized. OFFSET >= 0.

*A*

A is COMPLEX 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)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

**References:**

## subroutine claqps (integer M, integer N, integer OFFSET, integer NB, integer KB, complex, dimension( lda, * ) A, integer LDA, integer, dimension( * ) JPVT, complex, dimension( * ) TAU, real, dimension( * ) VN1, real, dimension( * ) VN2, complex, dimension( * ) AUXV, complex, dimension( ldf, * ) F, integer LDF)¶

**CLAQPS** computes a step of QR factorization with column
pivoting of a real m-by-n matrix A by using BLAS level 3.

**Purpose:**

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.

**Parameters**

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*OFFSET*

OFFSET is INTEGER

The number of rows of A that have been factorized in

previous steps.

*NB*

NB is INTEGER

The number of columns to factorize.

*KB*

KB is INTEGER

The number of columns actually factorized.

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Partial column norm updating strategy modified on April 2011 Z. Drmac and Z.
Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.

**References:**

## subroutine claqr0 (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, complex, dimension( ldh, * ) H, integer LDH, complex, dimension( * ) W, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, complex, dimension( * ) WORK, integer LWORK, integer INFO)¶

**CLAQR0** computes the eigenvalues of a Hessenberg matrix, and
optionally the matrices from the Schur decomposition.

**Purpose:**

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.

**Parameters**

*WANTT*

WANTT is LOGICAL

= .TRUE. : the full Schur form T is required;

= .FALSE.: only eigenvalues are required.

*WANTZ*

WANTZ is LOGICAL

= .TRUE. : the matrix of Schur vectors Z is required;

= .FALSE.: Schur vectors are not required.

*N*

N is INTEGER

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

*ILO*

ILO is INTEGER

*IHI*

IHI is INTEGER

It is assumed that H is already upper triangular in rows

and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,

H(ILO,ILO-1) is zero. ILO and IHI are normally set by a

previous call to 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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

**References:**

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR

Algorithm Part I: Maintaining Well Focused Shifts, and Level 3

Performance, SIAM Journal of Matrix Analysis, volume 23, pages

929--947, 2002.

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

## subroutine claqr1 (integer N, complex, dimension( ldh, * ) H, integer LDH, complex S1, complex S2, complex, dimension( * ) V)¶

**CLAQR1** sets a scalar multiple of the first column of the
product of 2-by-2 or 3-by-3 matrix H and specified shifts.

**Purpose:**

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

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine claqr2 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT, integer NW, complex, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, complex, dimension( * ) SH, complex, dimension( ldv, * ) V, integer LDV, integer NH, complex, dimension( ldt, * ) T, integer LDT, integer NV, complex, dimension( ldwv, * ) WV, integer LDWV, complex, dimension( * ) WORK, integer LWORK)¶

**CLAQR2** performs the unitary similarity transformation of a
Hessenberg matrix to detect and deflate fully converged eigenvalues from a
trailing principal submatrix (aggressive early deflation).

**Purpose:**

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.

**Parameters**

*WANTT*

WANTT is LOGICAL

If .TRUE., then the Hessenberg matrix H is fully updated

so that the triangular Schur factor may be

computed (in cooperation with the calling subroutine).

If .FALSE., then only enough of H is updated to preserve

the eigenvalues.

*WANTZ*

WANTZ is LOGICAL

If .TRUE., then the unitary matrix Z is updated so

so that the unitary Schur factor may be computed

(in cooperation with the calling subroutine).

If .FALSE., then Z is not referenced.

*N*

N is INTEGER

The order of the matrix H and (if WANTZ is .TRUE.) the

order of the unitary matrix Z.

*KTOP*

KTOP is INTEGER

It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.

KBOT and KTOP together determine an isolated block

along the diagonal of the Hessenberg matrix.

*KBOT*

KBOT is INTEGER

It is assumed without a check that either

KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together

determine an isolated block along the diagonal of the

Hessenberg matrix.

*NW*

NW is INTEGER

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

*H*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine claqr3 (logical WANTT, logical WANTZ, integer N, integer KTOP, integer KBOT, integer NW, complex, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, integer NS, integer ND, complex, dimension( * ) SH, complex, dimension( ldv, * ) V, integer LDV, integer NH, complex, dimension( ldt, * ) T, integer LDT, integer NV, complex, dimension( ldwv, * ) WV, integer LDWV, complex, dimension( * ) WORK, integer LWORK)¶

**CLAQR3** performs the unitary similarity transformation of a
Hessenberg matrix to detect and deflate fully converged eigenvalues from a
trailing principal submatrix (aggressive early deflation).

**Purpose:**

Aggressive early deflation:

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.

**Parameters**

*WANTT*

WANTT is LOGICAL

If .TRUE., then the Hessenberg matrix H is fully updated

so that the triangular Schur factor may be

computed (in cooperation with the calling subroutine).

If .FALSE., then only enough of H is updated to preserve

the eigenvalues.

*WANTZ*

WANTZ is LOGICAL

If .TRUE., then the unitary matrix Z is updated so

so that the unitary Schur factor may be computed

(in cooperation with the calling subroutine).

If .FALSE., then Z is not referenced.

*N*

N is INTEGER

The order of the matrix H and (if WANTZ is .TRUE.) the

order of the unitary matrix Z.

*KTOP*

KTOP is INTEGER

It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.

KBOT and KTOP together determine an isolated block

along the diagonal of the Hessenberg matrix.

*KBOT*

KBOT is INTEGER

It is assumed without a check that either

KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together

determine an isolated block along the diagonal of the

Hessenberg matrix.

*NW*

NW is INTEGER

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

*H*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## subroutine claqr4 (logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, complex, dimension( ldh, * ) H, integer LDH, complex, dimension( * ) W, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, complex, dimension( * ) WORK, integer LWORK, integer INFO)¶

**CLAQR4** computes the eigenvalues of a Hessenberg matrix, and
optionally the matrices from the Schur decomposition.

**Purpose:**

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.

**Parameters**

*WANTT*

WANTT is LOGICAL

= .TRUE. : the full Schur form T is required;

= .FALSE.: only eigenvalues are required.

*WANTZ*

WANTZ is LOGICAL

= .TRUE. : the matrix of Schur vectors Z is required;

= .FALSE.: Schur vectors are not required.

*N*

N is INTEGER

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

*ILO*

ILO is INTEGER

*IHI*

IHI is INTEGER

It is assumed that H is already upper triangular in rows

and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,

H(ILO,ILO-1) is zero. ILO and IHI are normally set by a

previous call to 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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

**References:**

K. Braman, R. Byers and R. Mathias, The Multi-Shift QR

Algorithm Part I: Maintaining Well Focused Shifts, and Level 3

Performance, SIAM Journal of Matrix Analysis, volume 23, pages

929--947, 2002.

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

## subroutine claqr5 (logical WANTT, logical WANTZ, integer KACC22, integer N, integer KTOP, integer KBOT, integer NSHFTS, complex, dimension( * ) S, complex, dimension( ldh, * ) H, integer LDH, integer ILOZ, integer IHIZ, complex, dimension( ldz, * ) Z, integer LDZ, complex, dimension( ldv, * ) V, integer LDV, complex, dimension( ldu, * ) U, integer LDU, integer NV, complex, dimension( ldwv, * ) WV, integer LDWV, integer NH, complex, dimension( ldwh, * ) WH, integer LDWH)¶

**CLAQR5** performs a single small-bulge multi-shift QR
sweep.

**Purpose:**

CLAQR5 called by CLAQR0 performs a

single small-bulge multi-shift QR sweep.

**Parameters**

*WANTT*

WANTT is LOGICAL

WANTT = .true. if the triangular Schur factor

is being computed. WANTT is set to .false. otherwise.

*WANTZ*

WANTZ is LOGICAL

WANTZ = .true. if the unitary Schur factor is being

computed. WANTZ is set to .false. otherwise.

*KACC22*

KACC22 is INTEGER with value 0, 1, or 2.

Specifies the computation mode of far-from-diagonal

orthogonal updates.

= 0: 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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Lars Karlsson, Daniel Kressner, and Bruno Lang

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

**References:**

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

## subroutine claqsb (character UPLO, integer N, integer KD, complex, dimension( ldab, * ) AB, integer LDAB, real, dimension( * ) S, real SCOND, real AMAX, character EQUED)¶

**CLAQSB** scales a symmetric/Hermitian band matrix, using
scaling factors computed by spbequ.

**Purpose:**

CLAQSB equilibrates a symmetric band matrix A using the scaling

factors in the vector S.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the

symmetric matrix A is stored.

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

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

*KD*

KD is INTEGER

The number of super-diagonals of the matrix A if UPLO = 'U',

or the number of sub-diagonals if UPLO = 'L'. KD >= 0.

*AB*

AB is COMPLEX 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).

**Internal Parameters:**

THRESH is a threshold value used to decide if scaling should be done

based on the ratio of the scaling factors. If SCOND < THRESH,

scaling is done.

LARGE and SMALL are threshold values used to decide if scaling should

be done based on the absolute size of the largest matrix element.

If AMAX > LARGE or AMAX < SMALL, scaling is done.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claqsp (character UPLO, integer N, complex, dimension( * ) AP, real, dimension( * ) S, real SCOND, real AMAX, character EQUED)¶

**CLAQSP** scales a symmetric/Hermitian matrix in packed
storage, using scaling factors computed by sppequ.

**Purpose:**

CLAQSP equilibrates a symmetric matrix A using the scaling factors

in the vector S.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the

symmetric matrix A is stored.

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

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

*AP*

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

**Internal Parameters:**

THRESH is a threshold value used to decide if scaling should be done

based on the ratio of the scaling factors. If SCOND < THRESH,

scaling is done.

LARGE and SMALL are threshold values used to decide if scaling should

be done based on the absolute size of the largest matrix element.

If AMAX > LARGE or AMAX < SMALL, scaling is done.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clar1v (integer N, integer B1, integer BN, real LAMBDA, real, dimension( * ) D, real, dimension( * ) L, real, dimension( * ) LD, real, dimension( * ) LLD, real PIVMIN, real GAPTOL, complex, dimension( * ) Z, logical WANTNC, integer NEGCNT, real ZTZ, real MINGMA, integer R, integer, dimension( * ) ISUPPZ, real NRMINV, real RESID, real RQCORR, real, dimension( * ) WORK)¶

**CLAR1V** computes the (scaled) r-th column of the inverse of
the submatrix in rows b1 through bn of the tridiagonal matrix LDLT -
λI.

**Purpose:**

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.

**Parameters**

*N*

N is INTEGER

The order of the matrix L D L**T.

*B1*

B1 is INTEGER

First index of the submatrix of L D L**T.

*BN*

BN is INTEGER

Last index of the submatrix of L D L**T.

*LAMBDA*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Jim Demmel, University of California, Berkeley, USA

Inderjit Dhillon, University of Texas, Austin, USA

Osni Marques, LBNL/NERSC, USA

Christof Voemel, University of California, Berkeley, USA

## subroutine clar2v (integer N, complex, dimension( * ) X, complex, dimension( * ) Y, complex, dimension( * ) Z, integer INCX, real, dimension( * ) C, complex, dimension( * ) S, integer INCC)¶

**CLAR2V** applies a vector of plane rotations with real
cosines and complex sines from both sides to a sequence of 2-by-2
symmetric/Hermitian matrices.

**Purpose:**

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

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clarcm (integer M, integer N, real, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldc, * ) C, integer LDC, real, dimension( * ) RWORK)¶

**CLARCM** copies all or part of a real two-dimensional array
to a complex array.

**Purpose:**

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.

**Parameters**

*M*

M is INTEGER

The number of rows of the matrix A and of the matrix C.

M >= 0.

*N*

N is INTEGER

The number of columns and rows of the matrix B and

the number of columns of the matrix C.

N >= 0.

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clarf (character SIDE, integer M, integer N, complex, dimension( * ) V, integer INCV, complex TAU, complex, dimension( ldc, * ) C, integer LDC, complex, dimension( * ) WORK)¶

**CLARF** applies an elementary reflector to a general
rectangular matrix.

**Purpose:**

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.

**Parameters**

*SIDE*

SIDE is CHARACTER*1

= 'L': form H * C

= 'R': form C * H

*M*

M is INTEGER

The number of rows of the matrix C.

*N*

N is INTEGER

The number of columns of the matrix C.

*V*

V is COMPLEX 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'

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clarfb (character SIDE, character TRANS, character DIRECT, character STOREV, integer M, integer N, integer K, complex, dimension( ldv, * ) V, integer LDV, complex, dimension( ldt, * ) T, integer LDT, complex, dimension( ldc, * ) C, integer LDC, complex, dimension( ldwork, * ) WORK, integer LDWORK)¶

**CLARFB** applies a block reflector or its conjugate-transpose
to a general rectangular matrix.

**Purpose:**

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.

**Parameters**

*SIDE*

SIDE is CHARACTER*1

= 'L': apply H or H**H from the Left

= 'R': apply H or H**H from the Right

*TRANS*

TRANS is CHARACTER*1

= 'N': apply H (No transpose)

= 'C': apply H**H (Conjugate transpose)

*DIRECT*

DIRECT is CHARACTER*1

Indicates how H is formed from a product of elementary

reflectors

= 'F': H = H(1) H(2) . . . H(k) (Forward)

= 'B': H = H(k) . . . H(2) H(1) (Backward)

*STOREV*

STOREV is CHARACTER*1

Indicates how the vectors which define the elementary

reflectors are stored:

= 'C': Columnwise

= 'R': Rowwise

*M*

M is INTEGER

The number of rows of the matrix C.

*N*

N is INTEGER

The number of columns of the matrix C.

*K*

K is INTEGER

The order of the matrix T (= the number of elementary

reflectors whose product defines the block reflector).

If SIDE = 'L', M >= K >= 0;

if SIDE = 'R', N >= K >= 0.

*V*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

The shape of the matrix V and the storage of the vectors which define

the H(i) is best illustrated by the following example with n = 5 and

k = 3. The elements equal to 1 are not stored; the corresponding

array elements are modified but restored on exit. The rest of the

array is not used.

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

V = ( 1 ) V = ( 1 v1 v1 v1 v1 )

( v1 1 ) ( 1 v2 v2 v2 )

( v1 v2 1 ) ( 1 v3 v3 )

( v1 v2 v3 )

( v1 v2 v3 )

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

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

( v1 v2 v3 ) ( v2 v2 v2 1 )

( 1 v2 v3 ) ( v3 v3 v3 v3 1 )

( 1 v3 )

( 1 )

## subroutine clarfb_gett (character IDENT, integer M, integer N, integer K, complex, dimension( ldt, * ) T, integer LDT, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldwork, * ) WORK, integer LDWORK)¶

**CLARFB_GETT**

**Purpose:**

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.

**Parameters**

*IDENT*

IDENT is CHARACTER*1

If IDENT = not 'I', or not 'i', then V1 is unit

lower-triangular and stored in the left K-by-K block of

the input matrix A,

If IDENT = 'I' or 'i', then V1 is an identity matrix and

not stored.

See Further Details section.

*M*

M is INTEGER

The number of rows of the matrix B.

M >= 0.

*N*

N is INTEGER

The number of columns of the matrices A and B.

N >= 0.

*K*

K is INTEGER

The number or rows of the matrix A.

K is also order of the matrix T, i.e. the number of

elementary reflectors whose product defines the block

reflector. 0 <= K <= N.

*T*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

November 2020, Igor Kozachenko,

Computer Science Division,

University of California, Berkeley

**Further Details:**

(1) Description of the Algebraic Operation.

The matrix A is a K-by-N matrix composed of two column block

matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):

A = ( A1, A2 ).

The matrix B is an M-by-N matrix composed of two column block

matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):

B = ( B1, B2 ).

Perform the operation:

( A_out ) := H * ( A_in ) = ( I - V * T * V**H ) * ( A_in ) =

( B_out ) ( B_in ) ( B_in )

= ( I - ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in )

( V2 ) ( B_in )

On input:

a) ( A_in ) consists of two block columns:

( B_in )

( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))

( B_in ) (( B1_in ) ( B2_in )) (( 0 ) ( B2_in )),

where the column blocks are:

( A1_in ) is a K-by-K upper-triangular matrix stored in the

upper triangular part of the array A(1:K,1:K).

( B1_in ) is an M-by-K rectangular ZERO matrix and not stored.

( A2_in ) is a K-by-(N-K) rectangular matrix stored

in the array A(1:K,K+1:N).

( B2_in ) is an M-by-(N-K) rectangular matrix stored

in the array B(1:M,K+1:N).

b) V = ( V1 )

( V2 )

where:

1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;

2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,

stored in the lower-triangular part of the array

A(1:K,1:K) (ones are not stored),

and V2 is an M-by-K rectangular stored the array B(1:M,1:K),

(because on input B1_in is a rectangular zero

matrix that is not stored and the space is

used to store V2).

c) T is a K-by-K upper-triangular matrix stored

in the array T(1:K,1:K).

On output:

a) ( A_out ) consists of two block columns:

( B_out )

( A_out ) = (( A1_out ) ( A2_out ))

( B_out ) (( B1_out ) ( B2_out )),

where the column blocks are:

( A1_out ) is a K-by-K square matrix, or a K-by-K

upper-triangular matrix, if V1 is an

identity matrix. AiOut is stored in

the array A(1:K,1:K).

( B1_out ) is an M-by-K rectangular matrix stored

in the array B(1:M,K:N).

( A2_out ) is a K-by-(N-K) rectangular matrix stored

in the array A(1:K,K+1:N).

( B2_out ) is an M-by-(N-K) rectangular matrix stored

in the array B(1:M,K+1:N).

The operation above can be represented as the same operation

on each block column:

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

( B1_out ) ( 0 ) ( 0 )

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

( B2_out ) ( B2_in ) ( B2_in )

If IDENT != 'I':

The computation for column block 1:

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

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

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

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

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

If IDENT == 'I':

The operation for column block 1:

A1_out: = A1_in - V1*T*A1_in

B1_out: = - V2*T*A1_in

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

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

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

(2) Description of the Algorithmic Computation.

In the first step, we compute column block 2, i.e. A2 and B2.

Here, we need to use the K-by-(N-K) rectangular workspace

matrix W2 that is of the same size as the matrix A2.

W2 is stored in the array WORK(1:K,1:(N-K)).

In the second step, we compute column block 1, i.e. A1 and B1.

Here, we need to use the K-by-K square workspace matrix W1

that is of the same size as the as the matrix A1.

W1 is stored in the array WORK(1:K,1:K).

NOTE: Hence, in this routine, we need the workspace array WORK

only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from

the first step and W1 from the second step.

Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I',

more computations than in the Case (B).

if( IDENT != 'I' ) then

if ( N > K ) then

(First Step - column block 2)

col2_(1) W2: = A2

col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2

col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2

col2_(4) W2: = T * W2

col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2

col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2

col2_(7) A2: = A2 - W2

else

(Second Step - column block 1)

col1_(1) W1: = A1

col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1

col1_(3) W1: = T * W1

col1_(4) B1: = - V2 * W1 = - B1 * W1

col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1

col1_(6) square A1: = A1 - W1

end if

end if

Case (B), when V1 is an identity matrix, i.e. IDENT == 'I',

less computations than in the Case (A)

if( IDENT == 'I' ) then

if ( N > K ) then

(First Step - column block 2)

col2_(1) W2: = A2

col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2

col2_(4) W2: = T * W2

col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2

col2_(7) A2: = A2 - W2

else

(Second Step - column block 1)

col1_(1) W1: = A1

col1_(3) W1: = T * W1

col1_(4) B1: = - V2 * W1 = - B1 * W1

col1_(6) upper-triangular_of_(A1): = A1 - W1

end if

end if

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

algorithm:

if ( N > K ) then

(First Step - column block 2)

col2_(1) W2: = A2

if( IDENT != 'I' ) then

col2_(2) W2: = (V1**H) * W2

= (unit_lower_tr_of_(A1)**H) * W2

end if

col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2]

col2_(4) W2: = T * W2

col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2

if( IDENT != 'I' ) then

col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2

end if

col2_(7) A2: = A2 - W2

else

(Second Step - column block 1)

col1_(1) W1: = A1

if( IDENT != 'I' ) then

col1_(2) W1: = (V1**H) * W1

= (unit_lower_tr_of_(A1)**H) * W1

end if

col1_(3) W1: = T * W1

col1_(4) B1: = - V2 * W1 = - B1 * W1

if( IDENT != 'I' ) then

col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1

col1_(6_a) below_diag_of_(A1): = - below_diag_of_(W1)

end if

col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)

end if

## subroutine clarfg (integer N, complex ALPHA, complex, dimension( * ) X, integer INCX, complex TAU)¶

**CLARFG** generates an elementary reflector (Householder
matrix).

**Purpose:**

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 .

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clarfgp (integer N, complex ALPHA, complex, dimension( * ) X, integer INCX, complex TAU)¶

**CLARFGP** generates an elementary reflector (Householder
matrix) with non-negative beta.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clarft (character DIRECT, character STOREV, integer N, integer K, complex, dimension( ldv, * ) V, integer LDV, complex, dimension( * ) TAU, complex, dimension( ldt, * ) T, integer LDT)¶

**CLARFT** forms the triangular factor T of a block reflector H
= I - vtvH

**Purpose:**

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

**Parameters**

*DIRECT*

DIRECT is CHARACTER*1

Specifies the order in which the elementary reflectors are

multiplied to form the block reflector:

= 'F': H = H(1) H(2) . . . H(k) (Forward)

= 'B': H = H(k) . . . H(2) H(1) (Backward)

*STOREV*

STOREV is CHARACTER*1

Specifies how the vectors which define the elementary

reflectors are stored (see also Further Details):

= 'C': columnwise

= 'R': rowwise

*N*

N is INTEGER

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

*K*

K is INTEGER

The order of the triangular factor T (= the number of

elementary reflectors). K >= 1.

*V*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

The shape of the matrix V and the storage of the vectors which define

the H(i) is best illustrated by the following example with n = 5 and

k = 3. The elements equal to 1 are not stored.

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

V = ( 1 ) V = ( 1 v1 v1 v1 v1 )

( v1 1 ) ( 1 v2 v2 v2 )

( v1 v2 1 ) ( 1 v3 v3 )

( v1 v2 v3 )

( v1 v2 v3 )

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

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

( v1 v2 v3 ) ( v2 v2 v2 1 )

( 1 v2 v3 ) ( v3 v3 v3 v3 1 )

( 1 v3 )

( 1 )

## subroutine clarfx (character SIDE, integer M, integer N, complex, dimension( * ) V, complex TAU, complex, dimension( ldc, * ) C, integer LDC, complex, dimension( * ) WORK)¶

**CLARFX** applies an elementary reflector to a general
rectangular matrix, with loop unrolling when the reflector has order
≤ 10.

**Purpose:**

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.

**Parameters**

*SIDE*

SIDE is CHARACTER*1

= 'L': form H * C

= 'R': form C * H

*M*

M is INTEGER

The number of rows of the matrix C.

*N*

N is INTEGER

The number of columns of the matrix C.

*V*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clarfy (character UPLO, integer N, complex, dimension( * ) V, integer INCV, complex TAU, complex, dimension( ldc, * ) C, integer LDC, complex, dimension( * ) WORK)¶

**CLARFY**

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the

Hermitian matrix C is stored.

= 'U': Upper triangle

= 'L': Lower triangle

*N*

N is INTEGER

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

*V*

V is COMPLEX 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)

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clargv (integer N, complex, dimension( * ) X, integer INCX, complex, dimension( * ) Y, integer INCY, real, dimension( * ) C, integer INCC)¶

**CLARGV** generates a vector of plane rotations with real
cosines and complex sines.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

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

This version has a few statements commented out for thread safety

(machine parameters are computed on each entry). 10 feb 03, SJH.

## subroutine clarnv (integer IDIST, integer, dimension( 4 ) ISEED, integer N, complex, dimension( * ) X)¶

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

**Purpose:**

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

normal distribution.

**Parameters**

*IDIST*

IDIST is INTEGER

Specifies the distribution of the random numbers:

= 1: real and imaginary parts each uniform (0,1)

= 2: real and imaginary parts each uniform (-1,1)

= 3: real and imaginary parts each normal (0,1)

= 4: uniformly distributed on the disc abs(z) < 1

= 5: uniformly distributed on the circle abs(z) = 1

*ISEED*

ISEED is INTEGER array, dimension (4)

On entry, the seed of the random number generator; the array

elements must be between 0 and 4095, and ISEED(4) must be

odd.

On exit, the seed is updated.

*N*

N is INTEGER

The number of random numbers to be generated.

*X*

X is COMPLEX array, dimension (N)

The generated random numbers.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

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.

## subroutine clarrv (integer N, real VL, real VU, real, dimension( * ) D, real, dimension( * ) L, real PIVMIN, integer, dimension( * ) ISPLIT, integer M, integer DOL, integer DOU, real MINRGP, real RTOL1, real RTOL2, real, dimension( * ) W, real, dimension( * ) WERR, real, dimension( * ) WGAP, integer, dimension( * ) IBLOCK, integer, dimension( * ) INDEXW, real, dimension( * ) GERS, complex, dimension( ldz, * ) Z, integer LDZ, integer, dimension( * ) ISUPPZ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO)¶

**CLARRV** computes the eigenvectors of the tridiagonal matrix
T = L D LT given L, D and the eigenvalues of L D LT.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

Jim Demmel, University of California, Berkeley, USA

Inderjit Dhillon, University of Texas, Austin, USA

Osni Marques, LBNL/NERSC, USA

Christof Voemel, University of California, Berkeley, USA

## subroutine clartv (integer N, complex, dimension( * ) X, integer INCX, complex, dimension( * ) Y, integer INCY, real, dimension( * ) C, complex, dimension( * ) S, integer INCC)¶

**CLARTV** applies a vector of plane rotations with real
cosines and complex sines to the elements of a pair of vectors.

**Purpose:**

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

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clascl (character TYPE, integer KL, integer KU, real CFROM, real CTO, integer M, integer N, complex, dimension( lda, * ) A, integer LDA, integer INFO)¶

**CLASCL** multiplies a general rectangular matrix by a real
scalar defined as cto/cfrom.

**Purpose:**

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.

**Parameters**

*TYPE*

TYPE is CHARACTER*1

TYPE indices the storage type of the input matrix.

= 'G': A is a full matrix.

= 'L': A is a lower triangular matrix.

= 'U': A is an upper triangular matrix.

= 'H': A is an upper Hessenberg matrix.

= 'B': A is a symmetric band matrix with lower bandwidth KL

and upper bandwidth KU and with the only the lower

half stored.

= 'Q': A is a symmetric band matrix with lower bandwidth KL

and upper bandwidth KU and with the only the upper

half stored.

= 'Z': A is a band matrix with lower bandwidth KL and upper

bandwidth KU. See 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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claset (character UPLO, integer M, integer N, complex ALPHA, complex BETA, complex, dimension( lda, * ) A, integer LDA)¶

**CLASET** initializes the off-diagonal elements and the
diagonal elements of a matrix to given values.

**Purpose:**

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

ALPHA on the offdiagonals.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies the part of the matrix A to be set.

= 'U': Upper triangular part is set. The lower triangle

is unchanged.

= 'L': Lower triangular part is set. The upper triangle

is unchanged.

Otherwise: All of the matrix A is set.

*M*

M is INTEGER

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

*N*

N is INTEGER

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

*ALPHA*

ALPHA is COMPLEX

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clasr (character SIDE, character PIVOT, character DIRECT, integer M, integer N, real, dimension( * ) C, real, dimension( * ) S, complex, dimension( lda, * ) A, integer LDA)¶

**CLASR** applies a sequence of plane rotations to a general
rectangular matrix.

**Purpose:**

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.

**Parameters**

*SIDE*

SIDE is CHARACTER*1

Specifies whether the plane rotation matrix P is applied to

A on the left or the right.

= 'L': Left, compute A := P*A

= 'R': Right, compute A:= A*P**T

*PIVOT*

PIVOT is CHARACTER*1

Specifies the plane for which P(k) is a plane rotation

matrix.

= 'V': Variable pivot, the plane (k,k+1)

= 'T': Top pivot, the plane (1,k+1)

= 'B': Bottom pivot, the plane (k,z)

*DIRECT*

DIRECT is CHARACTER*1

Specifies whether P is a forward or backward sequence of

plane rotations.

= 'F': Forward, P = P(z-1)*...*P(2)*P(1)

= 'B': Backward, P = P(1)*P(2)*...*P(z-1)

*M*

M is INTEGER

The number of rows of the matrix A. If m <= 1, an immediate

return is effected.

*N*

N is INTEGER

The number of columns of the matrix A. If n <= 1, an

immediate return is effected.

*C*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine claswp (integer N, complex, dimension( lda, * ) A, integer LDA, integer K1, integer K2, integer, dimension( * ) IPIV, integer INCX)¶

**CLASWP** performs a series of row interchanges on a general
rectangular matrix.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

Modified by

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

## subroutine clatbs (character UPLO, character TRANS, character DIAG, character NORMIN, integer N, integer KD, complex, dimension( ldab, * ) AB, integer LDAB, complex, dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)¶

**CLATBS** solves a triangular banded system of equations.

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the matrix A is upper or lower triangular.

= 'U': Upper triangular

= 'L': Lower triangular

*TRANS*

TRANS is CHARACTER*1

Specifies the operation applied to A.

= 'N': Solve A * x = s*b (No transpose)

= 'T': Solve A**T * x = s*b (Transpose)

= 'C': Solve A**H * x = s*b (Conjugate transpose)

*DIAG*

DIAG is CHARACTER*1

Specifies whether or not the matrix A is unit triangular.

= 'N': Non-unit triangular

= 'U': Unit triangular

*NORMIN*

NORMIN is CHARACTER*1

Specifies whether CNORM has been set or not.

= 'Y': CNORM contains the column norms on entry

= 'N': CNORM is not set on entry. On exit, the norms will

be computed and stored in CNORM.

*N*

N is INTEGER

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

*KD*

KD is INTEGER

The number of subdiagonals or superdiagonals in the

triangular matrix A. KD >= 0.

*AB*

AB is COMPLEX 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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

A rough bound on x is computed; if that is less than overflow, 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).

## subroutine clatdf (integer IJOB, integer N, complex, dimension( ldz, * ) Z, integer LDZ, complex, dimension( * ) RHS, real RDSUM, real RDSCAL, integer, dimension( * ) IPIV, integer, dimension( * ) JPIV)¶

**CLATDF** uses the LU factorization of the n-by-n matrix
computed by sgetc2 and computes a contribution to the reciprocal
Dif-estimate.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

**Contributors:**

**References:**

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

## subroutine clatps (character UPLO, character TRANS, character DIAG, character NORMIN, integer N, complex, dimension( * ) AP, complex, dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)¶

**CLATPS** solves a triangular system of equations with the
matrix held in packed storage.

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the matrix A is upper or lower triangular.

= 'U': Upper triangular

= 'L': Lower triangular

*TRANS*

TRANS is CHARACTER*1

Specifies the operation applied to A.

= 'N': Solve A * x = s*b (No transpose)

= 'T': Solve A**T * x = s*b (Transpose)

= 'C': Solve A**H * x = s*b (Conjugate transpose)

*DIAG*

DIAG is CHARACTER*1

Specifies whether or not the matrix A is unit triangular.

= 'N': Non-unit triangular

= 'U': Unit triangular

*NORMIN*

NORMIN is CHARACTER*1

Specifies whether CNORM has been set or not.

= 'Y': CNORM contains the column norms on entry

= 'N': CNORM is not set on entry. On exit, the norms will

be computed and stored in CNORM.

*N*

N is INTEGER

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

*AP*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

A rough bound on x is computed; if that is less than overflow, 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).

## subroutine clatrd (character UPLO, integer N, integer NB, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) E, complex, dimension( * ) TAU, complex, dimension( ldw, * ) W, integer LDW)¶

**CLATRD** reduces the first nb rows and columns of a
symmetric/Hermitian matrix A to real tridiagonal form by an unitary
similarity transformation.

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the upper or lower triangular part of the

Hermitian matrix A is stored:

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

The order of the matrix A.

*NB*

NB is INTEGER

The number of rows and columns to be reduced.

*A*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

If UPLO = 'U', the matrix Q is represented as a product of elementary

reflectors

Q = H(n) H(n-1) . . . H(n-nb+1).

Each H(i) has the form

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

where tau is a complex scalar, and v is a complex vector with

v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),

and tau in TAU(i-1).

If UPLO = 'L', the matrix Q is represented as a product of elementary

reflectors

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

Each H(i) has the form

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

where tau is a complex scalar, and v is a complex vector with

v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),

and tau in TAU(i).

The elements of the vectors v together form the n-by-nb matrix V

which is needed, with W, to apply the transformation to the unreduced

part of the matrix, using a Hermitian rank-2k update of the form:

A := A - V*W**H - W*V**H.

The contents of A on exit are illustrated by the following examples

with n = 5 and nb = 2:

if UPLO = 'U': if UPLO = 'L':

( a a a v4 v5 ) ( d )

( a a v4 v5 ) ( 1 d )

( a 1 v5 ) ( v1 1 a )

( d 1 ) ( v1 v2 a a )

( d ) ( v1 v2 a a a )

where d denotes a diagonal element of the reduced matrix, a denotes

an element of the original matrix that is unchanged, and vi denotes

an element of the vector defining H(i).

## subroutine clatrs (character UPLO, character TRANS, character DIAG, character NORMIN, integer N, complex, dimension( lda, * ) A, integer LDA, complex, dimension( * ) X, real SCALE, real, dimension( * ) CNORM, integer INFO)¶

**CLATRS** solves a triangular system of equations with the
scale factor set to prevent overflow.

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the matrix A is upper or lower triangular.

= 'U': Upper triangular

= 'L': Lower triangular

*TRANS*

TRANS is CHARACTER*1

Specifies the operation applied to A.

= 'N': Solve A * x = s*b (No transpose)

= 'T': Solve A**T * x = s*b (Transpose)

= 'C': Solve A**H * x = s*b (Conjugate transpose)

*DIAG*

DIAG is CHARACTER*1

Specifies whether or not the matrix A is unit triangular.

= 'N': Non-unit triangular

= 'U': Unit triangular

*NORMIN*

NORMIN is CHARACTER*1

Specifies whether CNORM has been set or not.

= 'Y': CNORM contains the column norms on entry

= 'N': CNORM is not set on entry. On exit, the norms will

be computed and stored in CNORM.

*N*

N is INTEGER

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

*A*

A is COMPLEX 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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

A rough bound on x is computed; if that is less than overflow, 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).

## subroutine clauu2 (character UPLO, integer N, complex, dimension( lda, * ) A, integer LDA, integer INFO)¶

**CLAUU2** computes the product UUH or LHL, where U and L are
upper or lower triangular matrices (unblocked algorithm).

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the triangular factor stored in the array A

is upper or lower triangular:

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

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

*A*

A is COMPLEX array, dimension (LDA,N)

On entry, the triangular factor U or L.

On exit, if UPLO = 'U', the upper triangle of A is

overwritten with the upper triangle of the product U * U**H;

if UPLO = 'L', the lower triangle of A is overwritten with

the lower triangle of the product L**H * L.

*LDA*

LDA is INTEGER

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

*INFO*

INFO is INTEGER

= 0: successful exit

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine clauum (character UPLO, integer N, complex, dimension( lda, * ) A, integer LDA, integer INFO)¶

**CLAUUM** computes the product UUH or LHL, where U and L are
upper or lower triangular matrices (blocked algorithm).

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

Specifies whether the triangular factor stored in the array A

is upper or lower triangular:

= 'U': Upper triangular

= 'L': Lower triangular

*N*

N is INTEGER

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

*A*

A is COMPLEX array, dimension (LDA,N)

On entry, the triangular factor U or L.

On exit, if UPLO = 'U', the upper triangle of A is

overwritten with the upper triangle of the product U * U**H;

if UPLO = 'L', the lower triangle of A is overwritten with

the lower triangle of the product L**H * L.

*LDA*

LDA is INTEGER

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

*INFO*

INFO is INTEGER

= 0: successful exit

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine crot (integer N, complex, dimension( * ) CX, integer INCX, complex, dimension( * ) CY, integer INCY, real C, complex S)¶

**CROT** applies a plane rotation with real cosine and complex
sine to a pair of complex vectors.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine cspmv (character UPLO, integer N, complex ALPHA, complex, dimension( * ) AP, complex, dimension( * ) X, integer INCX, complex BETA, complex, dimension( * ) Y, integer INCY)¶

**CSPMV** computes a matrix-vector product for complex vectors
using a complex symmetric packed matrix

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

On entry, UPLO specifies whether the upper or lower

triangular part of the matrix A is supplied in the packed

array AP as follows:

UPLO = 'U' or 'u' The upper triangular part of A is

supplied in AP.

UPLO = 'L' or 'l' The lower triangular part of A is

supplied in AP.

Unchanged on exit.

*N*

N is INTEGER

On entry, N specifies the order of the matrix A.

N must be at least zero.

Unchanged on exit.

*ALPHA*

ALPHA is COMPLEX

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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine cspr (character UPLO, integer N, complex ALPHA, complex, dimension( * ) X, integer INCX, complex, dimension( * ) AP)¶

**CSPR** performs the symmetrical rank-1 update of a complex
symmetric packed matrix.

**Purpose:**

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.

**Parameters**

*UPLO*

UPLO is CHARACTER*1

On entry, UPLO specifies whether the upper or lower

triangular part of the matrix A is supplied in the packed

array AP as follows:

UPLO = 'U' or 'u' The upper triangular part of A is

supplied in AP.

UPLO = 'L' or 'l' The lower triangular part of A is

supplied in AP.

Unchanged on exit.

*N*

N is INTEGER

On entry, N specifies the order of the matrix A.

N must be at least zero.

Unchanged on exit.

*ALPHA*

ALPHA is COMPLEX

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.

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine csrscl (integer N, real SA, complex, dimension( * ) SX, integer INCX)¶

**CSRSCL** multiplies a vector by the reciprocal of a real
scalar.

**Purpose:**

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.

**Parameters**

*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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## subroutine ctprfb (character SIDE, character TRANS, character DIRECT, character STOREV, integer M, integer N, integer K, integer L, complex, dimension( ldv, * ) V, integer LDV, complex, dimension( ldt, * ) T, integer LDT, complex, dimension( lda, * ) A, integer LDA, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldwork, * ) WORK, integer LDWORK)¶

**CTPRFB** applies a real or complex 'triangular-pentagonal'
blocked reflector to a real or complex matrix, which is composed of two
blocks.

**Purpose:**

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.

**Parameters**

*SIDE*

SIDE is CHARACTER*1

= 'L': apply H or H**H from the Left

= 'R': apply H or H**H from the Right

*TRANS*

TRANS is CHARACTER*1

= 'N': apply H (No transpose)

= 'C': apply H**H (Conjugate transpose)

*DIRECT*

DIRECT is CHARACTER*1

Indicates how H is formed from a product of elementary

reflectors

= 'F': H = H(1) H(2) . . . H(k) (Forward)

= 'B': H = H(k) . . . H(2) H(1) (Backward)

*STOREV*

STOREV is CHARACTER*1

Indicates how the vectors which define the elementary

reflectors are stored:

= 'C': Columns

= 'R': Rows

*M*

M is INTEGER

The number of rows of the matrix B.

M >= 0.

*N*

N is INTEGER

The number of columns of the matrix B.

N >= 0.

*K*

K is INTEGER

The order of the matrix T, i.e. the number of elementary

reflectors whose product defines the block reflector.

K >= 0.

*L*

L is INTEGER

The order of the trapezoidal part of V.

K >= L >= 0. See Further Details.

*V*

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Further Details:**

The matrix C is a composite matrix formed from blocks A and B.

The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,

and if SIDE = 'L', A is of size K-by-N.

If SIDE = 'R' and DIRECT = 'F', C = [A B].

If SIDE = 'L' and DIRECT = 'F', C = [A]

[B].

If SIDE = 'R' and DIRECT = 'B', C = [B A].

If SIDE = 'L' and DIRECT = 'B', C = [B]

[A].

The pentagonal matrix V is composed of a rectangular block V1 and a

trapezoidal block V2. The size of the trapezoidal block is determined by

the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;

if L=0, there is no trapezoidal block, thus V = V1 is rectangular.

If DIRECT = 'F' and STOREV = 'C': V = [V1]

[V2]

- V2 is upper trapezoidal (first L rows of K-by-K upper triangular)

If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]

- V2 is lower trapezoidal (first L columns of K-by-K lower triangular)

If DIRECT = 'B' and STOREV = 'C': V = [V2]

[V1]

- V2 is lower trapezoidal (last L rows of K-by-K lower triangular)

If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]

- V2 is upper trapezoidal (last L columns of K-by-K upper triangular)

If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.

If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.

If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.

If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.

## integer function icmax1 (integer N, complex, dimension(*) CX, integer INCX)¶

**ICMAX1** finds the index of the first vector element of
maximum absolute value.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## integer function ilaclc (integer M, integer N, complex, dimension( lda, * ) A, integer LDA)¶

**ILACLC** scans a matrix for its last non-zero column.

**Purpose:**

ILACLC scans A for its last non-zero column.

**Parameters**

*M*

M is INTEGER

The number of rows of the matrix A.

*N*

N is INTEGER

The number of columns of the matrix A.

*A*

A is COMPLEX array, dimension (LDA,N)

The m by n matrix A.

*LDA*

LDA is INTEGER

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## integer function ilaclr (integer M, integer N, complex, dimension( lda, * ) A, integer LDA)¶

**ILACLR** scans a matrix for its last non-zero row.

**Purpose:**

ILACLR scans A for its last non-zero row.

**Parameters**

*M*

M is INTEGER

The number of rows of the matrix A.

*N*

N is INTEGER

The number of columns of the matrix A.

*A*

A is COMPLEX array, dimension (LDA,N)

The m by n matrix A.

*LDA*

LDA is INTEGER

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

## integer function izmax1 (integer N, complex*16, dimension(*) ZX, integer INCX)¶

**IZMAX1** finds the index of the first vector element of
maximum absolute value.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

## real function scsum1 (integer N, complex, dimension( * ) CX, integer INCX)¶

**SCSUM1** forms the 1-norm of the complex vector using the
true absolute value.

**Purpose:**

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.

**Parameters**

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

**Author**

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

**Contributors:**

# Author¶

Generated automatically by Doxygen for LAPACK from the source code.

Wed Sep 1 2021 | Version 3.10.0 |