ZBBCSD computes the CS decomposition of a unitary matrix in


 bidiagonal-block form,


     [ B11 | B12 0  0 ]


     [  0  |  0 -I  0 ]


 X = [----------------]


     [ B21 | B22 0  0 ]


     [  0  |  0  0  I ]


                               [  C | -S  0  0 ]


                   [ U1 |    ] [  0 |  0 -I  0 ] [ V1 |    ]**H


                 = [---------] [---------------] [---------]   .


                   [    | U2 ] [  S |  C  0  0 ] [    | V2 ]


                               [  0 |  0  0  I ]


 X is M-by-M, its top-left block is P-by-Q, and Q must be no larger


 than P, M-P, or M-Q. (If Q is not the smallest index, then X must be


 transposed and/or permuted. This can be done in constant time using


 the TRANS and SIGNS options. See ZUNCSD for details.)


 The bidiagonal matrices B11, B12, B21, and B22 are represented


 implicitly by angles THETA(1:Q) and PHI(1:Q-1).


 The unitary matrices U1, U2, V1T, and V2T are input/output.


 The input matrices are pre- or post-multiplied by the appropriate


 singular vector matrices.

JOBU1

          JOBU1 is CHARACTER


          = 'Y':      U1 is updated;


          otherwise:  U1 is not updated.

JOBU2

          JOBU2 is CHARACTER


          = 'Y':      U2 is updated;


          otherwise:  U2 is not updated.

JOBV1T

          JOBV1T is CHARACTER


          = 'Y':      V1T is updated;


          otherwise:  V1T is not updated.

JOBV2T

          JOBV2T is CHARACTER


          = 'Y':      V2T is updated;


          otherwise:  V2T is not updated.

TRANS

          TRANS is CHARACTER


          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major


                      order;


          otherwise:  X, U1, U2, V1T, and V2T are stored in column-


                      major order.

M

          M is INTEGER


          The number of rows and columns in X, the unitary matrix in


          bidiagonal-block form.

P

          P is INTEGER


          The number of rows in the top-left block of X. 0 <= P <= M.

Q

          Q is INTEGER


          The number of columns in the top-left block of X.


          0 <= Q <= MIN(P,M-P,M-Q).

THETA

          THETA is DOUBLE PRECISION array, dimension (Q)


          On entry, the angles THETA(1),...,THETA(Q) that, along with


          PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block


          form. On exit, the angles whose cosines and sines define the


          diagonal blocks in the CS decomposition.

PHI

          PHI is DOUBLE PRECISION array, dimension (Q-1)


          The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,


          THETA(Q), define the matrix in bidiagonal-block form.

U1

          U1 is COMPLEX*16 array, dimension (LDU1,P)


          On entry, a P-by-P matrix. On exit, U1 is postmultiplied


          by the left singular vector matrix common to [ B11 ; 0 ] and


          [ B12 0 0 ; 0 -I 0 0 ].

LDU1

          LDU1 is INTEGER


          The leading dimension of the array U1, LDU1 >= MAX(1,P).

U2

          U2 is COMPLEX*16 array, dimension (LDU2,M-P)


          On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is


          postmultiplied by the left singular vector matrix common to


          [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].

LDU2

          LDU2 is INTEGER


          The leading dimension of the array U2, LDU2 >= MAX(1,M-P).

V1T

          V1T is COMPLEX*16 array, dimension (LDV1T,Q)


          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied


          by the conjugate transpose of the right singular vector


          matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].

LDV1T

          LDV1T is INTEGER


          The leading dimension of the array V1T, LDV1T >= MAX(1,Q).

V2T

          V2T is COMPLEX*16 array, dimension (LDV2T,M-Q)


          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is


          premultiplied by the conjugate transpose of the right


          singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and


          [ B22 0 0 ; 0 0 I ].

LDV2T

          LDV2T is INTEGER


          The leading dimension of the array V2T, LDV2T >= MAX(1,M-Q).

B11D

          B11D is DOUBLE PRECISION array, dimension (Q)


          When ZBBCSD converges, B11D contains the cosines of THETA(1),


          ..., THETA(Q). If ZBBCSD fails to converge, then B11D


          contains the diagonal of the partially reduced top-left


          block.

B11E

          B11E is DOUBLE PRECISION array, dimension (Q-1)


          When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails


          to converge, then B11E contains the superdiagonal of the


          partially reduced top-left block.

B12D

          B12D is DOUBLE PRECISION array, dimension (Q)


          When ZBBCSD converges, B12D contains the negative sines of


          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then


          B12D contains the diagonal of the partially reduced top-right


          block.

B12E

          B12E is DOUBLE PRECISION array, dimension (Q-1)


          When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails


          to converge, then B12E contains the subdiagonal of the


          partially reduced top-right block.

B21D

          B21D is DOUBLE PRECISION array, dimension (Q)


          When ZBBCSD converges, B21D contains the negative sines of


          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then


          B21D contains the diagonal of the partially reduced bottom-left


          block.

B21E

          B21E is DOUBLE PRECISION array, dimension (Q-1)


          When ZBBCSD converges, B21E contains zeros. If ZBBCSD fails


          to converge, then B21E contains the subdiagonal of the


          partially reduced bottom-left block.

B22D

          B22D is DOUBLE PRECISION array, dimension (Q)


          When ZBBCSD converges, B22D contains the negative sines of


          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then


          B22D contains the diagonal of the partially reduced bottom-right


          block.

B22E

          B22E is DOUBLE PRECISION array, dimension (Q-1)


          When ZBBCSD converges, B22E contains zeros. If ZBBCSD fails


          to converge, then B22E contains the subdiagonal of the


          partially reduced bottom-right block.

RWORK

          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))


          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.

LRWORK

          LRWORK is INTEGER


          The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).


          If LRWORK = -1, then a workspace query is assumed; the


          routine only calculates the optimal size of the RWORK array,


          returns this value as the first entry of the work array, and


          no error message related to LRWORK is issued by XERBLA.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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


          > 0:  if ZBBCSD did not converge, INFO specifies the number


                of nonzero entries in PHI, and B11D, B11E, etc.,


                contain the partially reduced matrix.

  TOLMUL  DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))


          TOLMUL controls the convergence criterion of the QR loop.


          Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they


          are within TOLMUL*EPS of either bound.

[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZBDSQR computes the singular values and, optionally, the right and/or


 left singular vectors from the singular value decomposition (SVD) of


 a real N-by-N (upper or lower) bidiagonal matrix B using the implicit


 zero-shift QR algorithm.  The SVD of B has the form


    B = Q * S * P**H


 where S is the diagonal matrix of singular values, Q is an orthogonal


 matrix of left singular vectors, and P is an orthogonal matrix of


 right singular vectors.  If left singular vectors are requested, this


 subroutine actually returns U*Q instead of Q, and, if right singular


 vectors are requested, this subroutine returns P**H*VT instead of


 P**H, for given complex input matrices U and VT.  When U and VT are


 the unitary matrices that reduce a general matrix A to bidiagonal


 form: A = U*B*VT, as computed by ZGEBRD, then


    A = (U*Q) * S * (P**H*VT)


 is the SVD of A.  Optionally, the subroutine may also compute Q**H*C


 for a given complex input matrix C.


 See 'Computing  Small Singular Values of Bidiagonal Matrices With


 Guaranteed High Relative Accuracy,' by J. Demmel and W. Kahan,


 LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,


 no. 5, pp. 873-912, Sept 1990) and


 'Accurate singular values and differential qd algorithms,' by


 B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics


 Department, University of California at Berkeley, July 1992


 for a detailed description of the algorithm.

UPLO

          UPLO is CHARACTER*1


          = 'U':  B is upper bidiagonal;


          = 'L':  B is lower bidiagonal.

N

          N is INTEGER


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

NCVT

          NCVT is INTEGER


          The number of columns of the matrix VT. NCVT >= 0.

NRU

          NRU is INTEGER


          The number of rows of the matrix U. NRU >= 0.

NCC

          NCC is INTEGER


          The number of columns of the matrix C. NCC >= 0.

D

          D is DOUBLE PRECISION array, dimension (N)


          On entry, the n diagonal elements of the bidiagonal matrix B.


          On exit, if INFO=0, the singular values of B in decreasing


          order.

E

          E is DOUBLE PRECISION array, dimension (N-1)


          On entry, the N-1 offdiagonal elements of the bidiagonal


          matrix B.


          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E


          will contain the diagonal and superdiagonal elements of a


          bidiagonal matrix orthogonally equivalent to the one given


          as input.

VT

          VT is COMPLEX*16 array, dimension (LDVT, NCVT)


          On entry, an N-by-NCVT matrix VT.


          On exit, VT is overwritten by P**H * VT.


          Not referenced if NCVT = 0.

LDVT

          LDVT is INTEGER


          The leading dimension of the array VT.


          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.

U

          U is COMPLEX*16 array, dimension (LDU, N)


          On entry, an NRU-by-N matrix U.


          On exit, U is overwritten by U * Q.


          Not referenced if NRU = 0.

LDU

          LDU is INTEGER


          The leading dimension of the array U.  LDU >= max(1,NRU).

C

          C is COMPLEX*16 array, dimension (LDC, NCC)


          On entry, an N-by-NCC matrix C.


          On exit, C is overwritten by Q**H * C.


          Not referenced if NCC = 0.

LDC

          LDC is INTEGER


          The leading dimension of the array C.


          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.

RWORK

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

INFO

          INFO is INTEGER


          = 0:  successful exit


          < 0:  If INFO = -i, the i-th argument had an illegal value


          > 0:  the algorithm did not converge; D and E contain the


                elements of a bidiagonal matrix which is orthogonally


                similar to the input matrix B;  if INFO = i, i


                elements of E have not converged to zero.

  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))


          TOLMUL controls the convergence criterion of the QR loop.


          If it is positive, TOLMUL*EPS is the desired relative


             precision in the computed singular values.


          If it is negative, abs(TOLMUL*EPS*sigma_max) is the


             desired absolute accuracy in the computed singular


             values (corresponds to relative accuracy


             abs(TOLMUL*EPS) in the largest singular value.


          abs(TOLMUL) should be between 1 and 1/EPS, and preferably


             between 10 (for fast convergence) and .1/EPS


             (for there to be some accuracy in the results).


          Default is to lose at either one eighth or 2 of the


             available decimal digits in each computed singular value


             (whichever is smaller).


  MAXITR  INTEGER, default = 6


          MAXITR controls the maximum number of passes of the


          algorithm through its inner loop. The algorithms stops


          (and so fails to converge) if the number of passes


          through the inner loop exceeds MAXITR*N**2.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZGGHD3 reduces a pair of complex matrices (A,B) to generalized upper


 Hessenberg form using unitary transformations, where A is a


 general matrix and B is upper triangular.  The form of the


 generalized eigenvalue problem is


    A*x = lambda*B*x,


 and B is typically made upper triangular by computing its QR


 factorization and moving the unitary matrix Q to the left side


 of the equation.


 This subroutine simultaneously reduces A to a Hessenberg matrix H:


    Q**H*A*Z = H


 and transforms B to another upper triangular matrix T:


    Q**H*B*Z = T


 in order to reduce the problem to its standard form


    H*y = lambda*T*y


 where y = Z**H*x.


 The unitary matrices Q and Z are determined as products of Givens


 rotations.  They may either be formed explicitly, or they may be


 postmultiplied into input matrices Q1 and Z1, so that


      Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H


      Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H


 If Q1 is the unitary matrix from the QR factorization of B in the


 original equation A*x = lambda*B*x, then ZGGHD3 reduces the original


 problem to generalized Hessenberg form.


 This is a blocked variant of CGGHRD, using matrix-matrix


 multiplications for parts of the computation to enhance performance.

COMPQ

          COMPQ is CHARACTER*1


          = 'N': do not compute Q;


          = 'I': Q is initialized to the unit matrix, and the


                 unitary matrix Q is returned;


          = 'V': Q must contain a unitary matrix Q1 on entry,


                 and the product Q1*Q is returned.

COMPZ

          COMPZ is CHARACTER*1


          = 'N': do not compute Z;


          = 'I': Z is initialized to the unit matrix, and the


                 unitary matrix Z is returned;


          = 'V': Z must contain a unitary matrix Z1 on entry,


                 and the product Z1*Z is returned.

N

          N is INTEGER


          The order of the matrices A and B.  N >= 0.

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


          ILO and IHI mark the rows and columns of A which are to be


          reduced.  It is assumed that A is already upper triangular


          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are


          normally set by a previous call to ZGGBAL; otherwise they


          should be set to 1 and N respectively.


          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

A

          A is COMPLEX*16 array, dimension (LDA, N)


          On entry, the N-by-N general matrix to be reduced.


          On exit, the upper triangle and the first subdiagonal of A


          are overwritten with the upper Hessenberg matrix H, and the


          rest is set to zero.

LDA

          LDA is INTEGER


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

B

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


          On entry, the N-by-N upper triangular matrix B.


          On exit, the upper triangular matrix T = Q**H B Z.  The


          elements below the diagonal are set to zero.

LDB

          LDB is INTEGER


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

Q

          Q is COMPLEX*16 array, dimension (LDQ, N)


          On entry, if COMPQ = 'V', the unitary matrix Q1, typically


          from the QR factorization of B.


          On exit, if COMPQ='I', the unitary matrix Q, and if


          COMPQ = 'V', the product Q1*Q.


          Not referenced if COMPQ='N'.

LDQ

          LDQ is INTEGER


          The leading dimension of the array Q.


          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.

Z

          Z is COMPLEX*16 array, dimension (LDZ, N)


          On entry, if COMPZ = 'V', the unitary matrix Z1.


          On exit, if COMPZ='I', the unitary matrix Z, and if


          COMPZ = 'V', the product Z1*Z.


          Not referenced if COMPZ='N'.

LDZ

          LDZ is INTEGER


          The leading dimension of the array Z.


          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.

WORK

          WORK is COMPLEX*16 array, dimension (LWORK)


          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

LWORK

          LWORK is INTEGER


          The length of the array WORK.  LWORK >= 1.


          For optimum performance LWORK >= 6*N*NB, where NB is the


          optimal blocksize.


          If LWORK = -1, then a workspace query is assumed; the routine


          only calculates the optimal size of the WORK array, returns


          this value as the first entry of the WORK array, and no error


          message related to LWORK is issued by XERBLA.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  This routine reduces A to Hessenberg form and maintains B in triangular form


  using a blocked variant of Moler and Stewart's original algorithm,


  as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti


  (BIT 2008).

 ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper


 Hessenberg form using unitary transformations, where A is a


 general matrix and B is upper triangular.  The form of the


 generalized eigenvalue problem is


    A*x = lambda*B*x,


 and B is typically made upper triangular by computing its QR


 factorization and moving the unitary matrix Q to the left side


 of the equation.


 This subroutine simultaneously reduces A to a Hessenberg matrix H:


    Q**H*A*Z = H


 and transforms B to another upper triangular matrix T:


    Q**H*B*Z = T


 in order to reduce the problem to its standard form


    H*y = lambda*T*y


 where y = Z**H*x.


 The unitary matrices Q and Z are determined as products of Givens


 rotations.  They may either be formed explicitly, or they may be


 postmultiplied into input matrices Q1 and Z1, so that


      Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H


      Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H


 If Q1 is the unitary matrix from the QR factorization of B in the


 original equation A*x = lambda*B*x, then ZGGHRD reduces the original


 problem to generalized Hessenberg form.

COMPQ

          COMPQ is CHARACTER*1


          = 'N': do not compute Q;


          = 'I': Q is initialized to the unit matrix, and the


                 unitary matrix Q is returned;


          = 'V': Q must contain a unitary matrix Q1 on entry,


                 and the product Q1*Q is returned.

COMPZ

          COMPZ is CHARACTER*1


          = 'N': do not compute Z;


          = 'I': Z is initialized to the unit matrix, and the


                 unitary matrix Z is returned;


          = 'V': Z must contain a unitary matrix Z1 on entry,


                 and the product Z1*Z is returned.

N

          N is INTEGER


          The order of the matrices A and B.  N >= 0.

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


          ILO and IHI mark the rows and columns of A which are to be


          reduced.  It is assumed that A is already upper triangular


          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are


          normally set by a previous call to ZGGBAL; otherwise they


          should be set to 1 and N respectively.


          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

A

          A is COMPLEX*16 array, dimension (LDA, N)


          On entry, the N-by-N general matrix to be reduced.


          On exit, the upper triangle and the first subdiagonal of A


          are overwritten with the upper Hessenberg matrix H, and the


          rest is set to zero.

LDA

          LDA is INTEGER


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

B

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


          On entry, the N-by-N upper triangular matrix B.


          On exit, the upper triangular matrix T = Q**H B Z.  The


          elements below the diagonal are set to zero.

LDB

          LDB is INTEGER


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

Q

          Q is COMPLEX*16 array, dimension (LDQ, N)


          On entry, if COMPQ = 'V', the unitary matrix Q1, typically


          from the QR factorization of B.


          On exit, if COMPQ='I', the unitary matrix Q, and if


          COMPQ = 'V', the product Q1*Q.


          Not referenced if COMPQ='N'.

LDQ

          LDQ is INTEGER


          The leading dimension of the array Q.


          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.

Z

          Z is COMPLEX*16 array, dimension (LDZ, N)


          On entry, if COMPZ = 'V', the unitary matrix Z1.


          On exit, if COMPZ='I', the unitary matrix Z, and if


          COMPZ = 'V', the product Z1*Z.


          Not referenced if COMPZ='N'.

LDZ

          LDZ is INTEGER


          The leading dimension of the array Z.


          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  This routine reduces A to Hessenberg and B to triangular form by


  an unblocked reduction, as described in _Matrix_Computations_,


  by Golub and van Loan (Johns Hopkins Press).

 ZGGQRF computes a generalized QR factorization of an N-by-M matrix A


 and an N-by-P matrix B:


             A = Q*R,        B = Q*T*Z,


 where Q is an N-by-N unitary matrix, Z is a P-by-P unitary matrix,


 and R and T assume one of the forms:


 if N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N,


                 (  0  ) N-M                         N   M-N


                    M


 where R11 is upper triangular, and


 if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N-P,


                  P-N  N                           ( T21 ) P


                                                      P


 where T12 or T21 is upper triangular.


 In particular, if B is square and nonsingular, the GQR factorization


 of A and B implicitly gives the QR factorization of inv(B)*A:


              inv(B)*A = Z**H * (inv(T)*R)


 where inv(B) denotes the inverse of the matrix B, and Z**H denotes the


 conjugate transpose of matrix Z.

N

          N is INTEGER


          The number of rows of the matrices A and B. N >= 0.

M

          M is INTEGER


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

P

          P is INTEGER


          The number of columns of the matrix B.  P >= 0.

A

          A is COMPLEX*16 array, dimension (LDA,M)


          On entry, the N-by-M matrix A.


          On exit, the elements on and above the diagonal of the array


          contain the min(N,M)-by-M upper trapezoidal matrix R (R is


          upper triangular if N >= M); the elements below the diagonal,


          with the array TAUA, represent the unitary matrix Q as a


          product of min(N,M) elementary reflectors (see Further


          Details).

LDA

          LDA is INTEGER


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

TAUA

          TAUA is COMPLEX*16 array, dimension (min(N,M))


          The scalar factors of the elementary reflectors which


          represent the unitary matrix Q (see Further Details).

B

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


          On entry, the N-by-P matrix B.


          On exit, if N <= P, the upper triangle of the subarray


          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;


          if N > P, the elements on and above the (N-P)-th subdiagonal


          contain the N-by-P upper trapezoidal matrix T; the remaining


          elements, with the array TAUB, represent the unitary


          matrix Z as a product of elementary reflectors (see Further


          Details).

LDB

          LDB is INTEGER


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

TAUB

          TAUB is COMPLEX*16 array, dimension (min(N,P))


          The scalar factors of the elementary reflectors which


          represent the unitary matrix Z (see Further Details).

WORK

          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))


          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= max(1,N,M,P).


          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),


          where NB1 is the optimal blocksize for the QR factorization


          of an N-by-M matrix, NB2 is the optimal blocksize for the


          RQ factorization of an N-by-P matrix, and NB3 is the optimal


          blocksize for a call of ZUNMQR.


          If LWORK = -1, then a workspace query is assumed; the routine


          only calculates the optimal size of the WORK array, returns


          this value as the first entry of the WORK array, and no error


          message related to LWORK is issued by XERBLA.

INFO

          INFO is INTEGER


           = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The matrix Q is represented as a product of elementary reflectors


     Q = H(1) H(2) . . . H(k), where k = min(n,m).


  Each H(i) has the form


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


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


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


  and taua in TAUA(i).


  To form Q explicitly, use LAPACK subroutine ZUNGQR.


  To use Q to update another matrix, use LAPACK subroutine ZUNMQR.


  The matrix Z is represented as a product of elementary reflectors


     Z = H(1) H(2) . . . H(k), where k = min(n,p).


  Each H(i) has the form


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


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


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


  B(n-k+i,1:p-k+i-1), and taub in TAUB(i).


  To form Z explicitly, use LAPACK subroutine ZUNGRQ.


  To use Z to update another matrix, use LAPACK subroutine ZUNMRQ.

 ZGGRQF computes a generalized RQ factorization of an M-by-N matrix A


 and a P-by-N matrix B:


             A = R*Q,        B = Z*T*Q,


 where Q is an N-by-N unitary matrix, Z is a P-by-P unitary


 matrix, and R and T assume one of the forms:


 if M <= N,  R = ( 0  R12 ) M,   or if M > N,  R = ( R11 ) M-N,


                  N-M  M                           ( R21 ) N


                                                      N


 where R12 or R21 is upper triangular, and


 if P >= N,  T = ( T11 ) N  ,   or if P < N,  T = ( T11  T12 ) P,


                 (  0  ) P-N                         P   N-P


                    N


 where T11 is upper triangular.


 In particular, if B is square and nonsingular, the GRQ factorization


 of A and B implicitly gives the RQ factorization of A*inv(B):


              A*inv(B) = (R*inv(T))*Z**H


 where inv(B) denotes the inverse of the matrix B, and Z**H denotes the


 conjugate transpose of the matrix Z.

M

          M is INTEGER


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

P

          P is INTEGER


          The number of rows of the matrix B.  P >= 0.

N

          N is INTEGER


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

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, the M-by-N matrix A.


          On exit, if M <= N, the upper triangle of the subarray


          A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R;


          if M > N, the elements on and above the (M-N)-th subdiagonal


          contain the M-by-N upper trapezoidal matrix R; the remaining


          elements, with the array TAUA, 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,M).

TAUA

          TAUA is COMPLEX*16 array, dimension (min(M,N))


          The scalar factors of the elementary reflectors which


          represent the unitary matrix Q (see Further Details).

B

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


          On entry, the P-by-N matrix B.


          On exit, the elements on and above the diagonal of the array


          contain the min(P,N)-by-N upper trapezoidal matrix T (T is


          upper triangular if P >= N); the elements below the diagonal,


          with the array TAUB, represent the unitary matrix Z as a


          product of elementary reflectors (see Further Details).

LDB

          LDB is INTEGER


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

TAUB

          TAUB is COMPLEX*16 array, dimension (min(P,N))


          The scalar factors of the elementary reflectors which


          represent the unitary matrix Z (see Further Details).

WORK

          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))


          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= max(1,N,M,P).


          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),


          where NB1 is the optimal blocksize for the RQ factorization


          of an M-by-N matrix, NB2 is the optimal blocksize for the


          QR factorization of a P-by-N matrix, and NB3 is the optimal


          blocksize for a call of ZUNMRQ.


          If LWORK = -1, then a workspace query is assumed; the routine


          only calculates the optimal size of the WORK array, returns


          this value as the first entry of the WORK array, and no error


          message related to LWORK is issued by XERBLA.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The matrix Q is represented as a product of elementary reflectors


     Q = H(1) H(2) . . . H(k), where k = min(m,n).


  Each H(i) has the form


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


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


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


  A(m-k+i,1:n-k+i-1), and taua in TAUA(i).


  To form Q explicitly, use LAPACK subroutine ZUNGRQ.


  To use Q to update another matrix, use LAPACK subroutine ZUNMRQ.


  The matrix Z is represented as a product of elementary reflectors


     Z = H(1) H(2) . . . H(k), where k = min(p,n).


  Each H(i) has the form


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


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


  v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i),


  and taub in TAUB(i).


  To form Z explicitly, use LAPACK subroutine ZUNGQR.


  To use Z to update another matrix, use LAPACK subroutine ZUNMQR.

 ZGGSVP3 computes unitary matrices U, V and Q such that


                    N-K-L  K    L


  U**H*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;


                 L ( 0     0   A23 )


             M-K-L ( 0     0    0  )


                  N-K-L  K    L


         =     K ( 0    A12  A13 )  if M-K-L < 0;


             M-K ( 0     0   A23 )


                  N-K-L  K    L


  V**H*B*Q =   L ( 0     0   B13 )


             P-L ( 0     0    0  )


 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular


 upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,


 otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective


 numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H.


 This decomposition is the preprocessing step for computing the


 Generalized Singular Value Decomposition (GSVD), see subroutine


 ZGGSVD3.

JOBU

          JOBU is CHARACTER*1


          = 'U':  Unitary matrix U is computed;


          = 'N':  U is not computed.

JOBV

          JOBV is CHARACTER*1


          = 'V':  Unitary matrix V is computed;


          = 'N':  V is not computed.

JOBQ

          JOBQ is CHARACTER*1


          = 'Q':  Unitary matrix Q is computed;


          = 'N':  Q is not computed.

M

          M is INTEGER


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

P

          P is INTEGER


          The number of rows of the matrix B.  P >= 0.

N

          N is INTEGER


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

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, the M-by-N matrix A.


          On exit, A contains the triangular (or trapezoidal) matrix


          described in the Purpose section.

LDA

          LDA is INTEGER


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

B

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


          On entry, the P-by-N matrix B.


          On exit, B contains the triangular matrix described in


          the Purpose section.

LDB

          LDB is INTEGER


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

TOLA

          TOLA is DOUBLE PRECISION

TOLB

          TOLB is DOUBLE PRECISION


          TOLA and TOLB are the thresholds to determine the effective


          numerical rank of matrix B and a subblock of A. Generally,


          they are set to


             TOLA = MAX(M,N)*norm(A)*MAZHEPS,


             TOLB = MAX(P,N)*norm(B)*MAZHEPS.


          The size of TOLA and TOLB may affect the size of backward


          errors of the decomposition.

K

          K is INTEGER

L

          L is INTEGER


          On exit, K and L specify the dimension of the subblocks


          described in Purpose section.


          K + L = effective numerical rank of (A**H,B**H)**H.

U

          U is COMPLEX*16 array, dimension (LDU,M)


          If JOBU = 'U', U contains the unitary matrix U.


          If JOBU = 'N', U is not referenced.

LDU

          LDU is INTEGER


          The leading dimension of the array U. LDU >= max(1,M) if


          JOBU = 'U'; LDU >= 1 otherwise.

V

          V is COMPLEX*16 array, dimension (LDV,P)


          If JOBV = 'V', V contains the unitary matrix V.


          If JOBV = 'N', V is not referenced.

LDV

          LDV is INTEGER


          The leading dimension of the array V. LDV >= max(1,P) if


          JOBV = 'V'; LDV >= 1 otherwise.

Q

          Q is COMPLEX*16 array, dimension (LDQ,N)


          If JOBQ = 'Q', Q contains the unitary matrix Q.


          If JOBQ = 'N', Q is not referenced.

LDQ

          LDQ is INTEGER


          The leading dimension of the array Q. LDQ >= max(1,N) if


          JOBQ = 'Q'; LDQ >= 1 otherwise.

IWORK

          IWORK is INTEGER array, dimension (N)

RWORK

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

TAU

          TAU is COMPLEX*16 array, dimension (N)

WORK

          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))


          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If LWORK = -1, then a workspace query is assumed; the routine


          only calculates the optimal size of the WORK array, returns


          this value as the first entry of the WORK array, and no error


          message related to LWORK is issued by XERBLA.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The subroutine uses LAPACK subroutine ZGEQP3 for the QR factorization


  with column pivoting to detect the effective numerical rank of the


  a matrix. It may be replaced by a better rank determination strategy.


  ZGGSVP3 replaces the deprecated subroutine ZGGSVP.

 ZGSVJ0 is called from ZGESVJ as a pre-processor and that is its main


 purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but


 it does not check convergence (stopping criterion). Few tuning


 parameters (marked by [TP]) are available for the implementer.

JOBV

          JOBV is CHARACTER*1


          Specifies whether the output from this procedure is used


          to compute the matrix V:


          = 'V': the product of the Jacobi rotations is accumulated


                 by postmulyiplying the N-by-N array V.


                (See the description of V.)


          = 'A': the product of the Jacobi rotations is accumulated


                 by postmulyiplying the MV-by-N array V.


                (See the descriptions of MV and V.)


          = 'N': the Jacobi rotations are not accumulated.

M

          M is INTEGER


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

N

          N is INTEGER


          The number of columns of the input matrix A.


          M >= N >= 0.

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, M-by-N matrix A, such that A*diag(D) represents


          the input matrix.


          On exit,


          A_onexit * diag(D_onexit) represents the input matrix A*diag(D)


          post-multiplied by a sequence of Jacobi rotations, where the


          rotation threshold and the total number of sweeps are given in


          TOL and NSWEEP, respectively.


          (See the descriptions of D, TOL and NSWEEP.)

LDA

          LDA is INTEGER


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

D

          D is COMPLEX*16 array, dimension (N)


          The array D accumulates the scaling factors from the complex scaled


          Jacobi rotations.


          On entry, A*diag(D) represents the input matrix.


          On exit, A_onexit*diag(D_onexit) represents the input matrix


          post-multiplied by a sequence of Jacobi rotations, where the


          rotation threshold and the total number of sweeps are given in


          TOL and NSWEEP, respectively.


          (See the descriptions of A, TOL and NSWEEP.)

SVA

          SVA is DOUBLE PRECISION array, dimension (N)


          On entry, SVA contains the Euclidean norms of the columns of


          the matrix A*diag(D).


          On exit, SVA contains the Euclidean norms of the columns of


          the matrix A_onexit*diag(D_onexit).

MV

          MV is INTEGER


          If JOBV = 'A', then MV rows of V are post-multipled by a


                           sequence of Jacobi rotations.


          If JOBV = 'N',   then MV is not referenced.

V

          V is COMPLEX*16 array, dimension (LDV,N)


          If JOBV = 'V' then N rows of V are post-multipled by a


                           sequence of Jacobi rotations.


          If JOBV = 'A' then MV rows of V are post-multipled by a


                           sequence of Jacobi rotations.


          If JOBV = 'N',   then V is not referenced.

LDV

          LDV is INTEGER


          The leading dimension of the array V,  LDV >= 1.


          If JOBV = 'V', LDV >= N.


          If JOBV = 'A', LDV >= MV.

EPS

          EPS is DOUBLE PRECISION


          EPS = DLAMCH('Epsilon')

SFMIN

          SFMIN is DOUBLE PRECISION


          SFMIN = DLAMCH('Safe Minimum')

TOL

          TOL is DOUBLE PRECISION


          TOL is the threshold for Jacobi rotations. For a pair


          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is


          applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.

NSWEEP

          NSWEEP is INTEGER


          NSWEEP is the number of sweeps of Jacobi rotations to be


          performed.

WORK

          WORK is COMPLEX*16 array, dimension (LWORK)

LWORK

          LWORK is INTEGER


          LWORK is the dimension of WORK. LWORK >= M.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

ZGSVJ0 is used just to enable ZGESVJ to call a simplified version of itself to work on a submatrix of the original matrix.

Please report all bugs and send interesting test examples and comments to drmac@math.hr. Thank you.

 ZGSVJ1 is called from ZGESVJ as a pre-processor and that is its main


 purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but


 it targets only particular pivots and it does not check convergence


 (stopping criterion). Few tuning parameters (marked by [TP]) are


 available for the implementer.


 Further Details


 ~~~~~~~~~~~~~~~


 ZGSVJ1 applies few sweeps of Jacobi rotations in the column space of


 the input M-by-N matrix A. The pivot pairs are taken from the (1,2)


 off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The


 block-entries (tiles) of the (1,2) off-diagonal block are marked by the


 [x]'s in the following scheme:


    | *  *  * [x] [x] [x]|


    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.


    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.


    |[x] [x] [x] *  *  * |


    |[x] [x] [x] *  *  * |


    |[x] [x] [x] *  *  * |


 In terms of the columns of A, the first N1 columns are rotated 'against'


 the remaining N-N1 columns, trying to increase the angle between the


 corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is


 tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter.


 The number of sweeps is given in NSWEEP and the orthogonality threshold


 is given in TOL.

JOBV

          JOBV is CHARACTER*1


          Specifies whether the output from this procedure is used


          to compute the matrix V:


          = 'V': the product of the Jacobi rotations is accumulated


                 by postmulyiplying the N-by-N array V.


                (See the description of V.)


          = 'A': the product of the Jacobi rotations is accumulated


                 by postmulyiplying the MV-by-N array V.


                (See the descriptions of MV and V.)


          = 'N': the Jacobi rotations are not accumulated.

M

          M is INTEGER


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

N

          N is INTEGER


          The number of columns of the input matrix A.


          M >= N >= 0.

N1

          N1 is INTEGER


          N1 specifies the 2 x 2 block partition, the first N1 columns are


          rotated 'against' the remaining N-N1 columns of A.

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, M-by-N matrix A, such that A*diag(D) represents


          the input matrix.


          On exit,


          A_onexit * D_onexit represents the input matrix A*diag(D)


          post-multiplied by a sequence of Jacobi rotations, where the


          rotation threshold and the total number of sweeps are given in


          TOL and NSWEEP, respectively.


          (See the descriptions of N1, D, TOL and NSWEEP.)

LDA

          LDA is INTEGER


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

D

          D is COMPLEX*16 array, dimension (N)


          The array D accumulates the scaling factors from the fast scaled


          Jacobi rotations.


          On entry, A*diag(D) represents the input matrix.


          On exit, A_onexit*diag(D_onexit) represents the input matrix


          post-multiplied by a sequence of Jacobi rotations, where the


          rotation threshold and the total number of sweeps are given in


          TOL and NSWEEP, respectively.


          (See the descriptions of N1, A, TOL and NSWEEP.)

SVA

          SVA is DOUBLE PRECISION array, dimension (N)


          On entry, SVA contains the Euclidean norms of the columns of


          the matrix A*diag(D).


          On exit, SVA contains the Euclidean norms of the columns of


          the matrix onexit*diag(D_onexit).

MV

          MV is INTEGER


          If JOBV = 'A', then MV rows of V are post-multipled by a


                           sequence of Jacobi rotations.


          If JOBV = 'N',   then MV is not referenced.

V

          V is COMPLEX*16 array, dimension (LDV,N)


          If JOBV = 'V' then N rows of V are post-multipled by a


                           sequence of Jacobi rotations.


          If JOBV = 'A' then MV rows of V are post-multipled by a


                           sequence of Jacobi rotations.


          If JOBV = 'N',   then V is not referenced.

LDV

          LDV is INTEGER


          The leading dimension of the array V,  LDV >= 1.


          If JOBV = 'V', LDV >= N.


          If JOBV = 'A', LDV >= MV.

EPS

          EPS is DOUBLE PRECISION


          EPS = DLAMCH('Epsilon')

SFMIN

          SFMIN is DOUBLE PRECISION


          SFMIN = DLAMCH('Safe Minimum')

TOL

          TOL is DOUBLE PRECISION


          TOL is the threshold for Jacobi rotations. For a pair


          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is


          applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.

NSWEEP

          NSWEEP is INTEGER


          NSWEEP is the number of sweeps of Jacobi rotations to be


          performed.

WORK

          WORK is COMPLEX*16 array, dimension (LWORK)

LWORK

          LWORK is INTEGER


          LWORK is the dimension of WORK. LWORK >= M.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Zlatko Drmac (Zagreb, Croatia)

 ZHBGST reduces a complex Hermitian-definite banded generalized


 eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,


 such that C has the same bandwidth as A.


 B must have been previously factorized as S**H*S by ZPBSTF, using a


 split Cholesky factorization. A is overwritten by C = X**H*A*X, where


 X = S**(-1)*Q and Q is a unitary matrix chosen to preserve the


 bandwidth of A.

VECT

          VECT is CHARACTER*1


          = 'N':  do not form the transformation matrix X;


          = 'V':  form X.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


          The order of the matrices A and B.  N >= 0.

KA

          KA is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.

KB

          KB is INTEGER


          The number of superdiagonals of the matrix B if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          On entry, the upper or lower triangle of the Hermitian band


          matrix A, stored in the first ka+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(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;


          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).


          On exit, the transformed matrix X**H*A*X, stored in the same


          format as A.

LDAB

          LDAB is INTEGER


          The leading dimension of the array AB.  LDAB >= KA+1.

BB

          BB is COMPLEX*16 array, dimension (LDBB,N)


          The banded factor S from the split Cholesky factorization of


          B, as returned by ZPBSTF, stored in the first kb+1 rows of


          the array.

LDBB

          LDBB is INTEGER


          The leading dimension of the array BB.  LDBB >= KB+1.

X

          X is COMPLEX*16 array, dimension (LDX,N)


          If VECT = 'V', the n-by-n matrix X.


          If VECT = 'N', the array X is not referenced.

LDX

          LDX is INTEGER


          The leading dimension of the array X.


          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.

WORK

          WORK is COMPLEX*16 array, dimension (N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHBTRD reduces a complex Hermitian band matrix A to real symmetric


 tridiagonal form T by a unitary similarity transformation:


 Q**H * A * Q = T.

VECT

          VECT is CHARACTER*1


          = 'N':  do not form Q;


          = 'V':  form Q;


          = 'U':  update a matrix X, by forming X*Q.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          On entry, the upper or lower triangle of the Hermitian 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, the diagonal elements of AB are overwritten by the


          diagonal elements of the tridiagonal matrix T; if KD > 0, the


          elements on the first superdiagonal (if UPLO = 'U') or the


          first subdiagonal (if UPLO = 'L') are overwritten by the


          off-diagonal elements of T; the rest of AB is overwritten by


          values generated during the reduction.

LDAB

          LDAB is INTEGER


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

D

          D is DOUBLE PRECISION array, dimension (N)


          The diagonal elements of the tridiagonal matrix T.

E

          E is DOUBLE PRECISION array, dimension (N-1)


          The off-diagonal elements of the tridiagonal matrix T:


          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.

Q

          Q is COMPLEX*16 array, dimension (LDQ,N)


          On entry, if VECT = 'U', then Q must contain an N-by-N


          matrix X; if VECT = 'N' or 'V', then Q need not be set.


          On exit:


          if VECT = 'V', Q contains the N-by-N unitary matrix Q;


          if VECT = 'U', Q contains the product X*Q;


          if VECT = 'N', the array Q is not referenced.

LDQ

          LDQ is INTEGER


          The leading dimension of the array Q.


          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.

WORK

          WORK is COMPLEX*16 array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Modified by Linda Kaufman, Bell Labs.

 ZHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric


 tridiagonal form T by a unitary similarity transformation:


 Q**H * A * Q = T.

STAGE1

          STAGE1 is CHARACTER*1


          = 'N':  'No': to mention that the stage 1 of the reduction  


                  from dense to band using the zhetrd_he2hb routine


                  was not called before this routine to reproduce AB. 


                  In other term this routine is called as standalone. 


          = 'Y':  'Yes': to mention that the stage 1 of the 


                  reduction from dense to band using the zhetrd_he2hb 


                  routine has been called to produce AB (e.g., AB is


                  the output of zhetrd_he2hb.

VECT

          VECT is CHARACTER*1


          = 'N':  No need for the Housholder representation, 


                  and thus LHOUS is of size max(1, 4*N);


          = 'V':  the Householder representation is needed to 


                  either generate or to apply Q later on, 


                  then LHOUS is to be queried and computed.


                  (NOT AVAILABLE IN THIS RELEASE).

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          On entry, the upper or lower triangle of the Hermitian 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, the diagonal elements of AB are overwritten by the


          diagonal elements of the tridiagonal matrix T; if KD > 0, the


          elements on the first superdiagonal (if UPLO = 'U') or the


          first subdiagonal (if UPLO = 'L') are overwritten by the


          off-diagonal elements of T; the rest of AB is overwritten by


          values generated during the reduction.

LDAB

          LDAB is INTEGER


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

D

          D is DOUBLE PRECISION array, dimension (N)


          The diagonal elements of the tridiagonal matrix T.

E

          E is DOUBLE PRECISION array, dimension (N-1)


          The off-diagonal elements of the tridiagonal matrix T:


          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.

HOUS

          HOUS is COMPLEX*16 array, dimension LHOUS, that


          store the Householder representation.

LHOUS

          LHOUS is INTEGER


          The dimension of the array HOUS. LHOUS = MAX(1, dimension)


          If LWORK = -1, or LHOUS=-1,


          then a query is assumed; the routine


          only calculates the optimal size of the HOUS array, returns


          this value as the first entry of the HOUS array, and no error


          message related to LHOUS is issued by XERBLA.


          LHOUS = MAX(1, dimension) where


          dimension = 4*N if VECT='N'


          not available now if VECT='H'     

WORK

          WORK is COMPLEX*16 array, dimension LWORK.

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK = MAX(1, dimension)


          If LWORK = -1, or LHOUS=-1,


          then a workspace query is assumed; the routine


          only calculates the optimal size of the WORK array, returns


          this value as the first entry of the WORK array, and no error


          message related to LWORK is issued by XERBLA.


          LWORK = MAX(1, dimension) where


          dimension   = (2KD+1)*N + KD*NTHREADS


          where KD is the blocking size of the reduction,


          FACTOPTNB is the blocking used by the QR or LQ


          algorithm, usually FACTOPTNB=128 is a good choice


          NTHREADS is the number of threads used when


          openMP compilation is enabled, otherwise =1.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Implemented by Azzam Haidar.


  All details are available on technical report, SC11, SC13 papers.


  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.


  Parallel reduction to condensed forms for symmetric eigenvalue problems


  using aggregated fine-grained and memory-aware kernels. In Proceedings


  of 2011 International Conference for High Performance Computing,


  Networking, Storage and Analysis (SC '11), New York, NY, USA,


  Article 8 , 11 pages.


  http://doi.acm.org/10.1145/2063384.2063394


  A. Haidar, J. Kurzak, P. Luszczek, 2013.


  An improved parallel singular value algorithm and its implementation 


  for multicore hardware, In Proceedings of 2013 International Conference


  for High Performance Computing, Networking, Storage and Analysis (SC '13).


  Denver, Colorado, USA, 2013.


  Article 90, 12 pages.


  http://doi.acm.org/10.1145/2503210.2503292


  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.


  A novel hybrid CPU-GPU generalized eigensolver for electronic structure 


  calculations based on fine-grained memory aware tasks.


  International Journal of High Performance Computing Applications.


  Volume 28 Issue 2, Pages 196-209, May 2014.


  http://hpc.sagepub.com/content/28/2/196 

 Level 3 BLAS like routine for C in RFP Format.


 ZHFRK performs one of the Hermitian rank--k operations


    C := alpha*A*A**H + beta*C,


 or


    C := alpha*A**H*A + beta*C,


 where alpha and beta are real scalars, C is an n--by--n Hermitian


 matrix and A is an n--by--k matrix in the first case and a k--by--n


 matrix in the second case.

TRANSR

          TRANSR is CHARACTER*1


          = 'N':  The Normal Form of RFP A is stored;


          = 'C':  The Conjugate-transpose Form of RFP A is stored.

UPLO

          UPLO is CHARACTER*1


           On  entry,   UPLO  specifies  whether  the  upper  or  lower


           triangular  part  of the  array  C  is to be  referenced  as


           follows:


              UPLO = 'U' or 'u'   Only the  upper triangular part of  C


                                  is to be referenced.


              UPLO = 'L' or 'l'   Only the  lower triangular part of  C


                                  is to be referenced.


           Unchanged on exit.

TRANS

          TRANS is CHARACTER*1


           On entry,  TRANS  specifies the operation to be performed as


           follows:


              TRANS = 'N' or 'n'   C := alpha*A*A**H + beta*C.


              TRANS = 'C' or 'c'   C := alpha*A**H*A + beta*C.


           Unchanged on exit.

N

          N is INTEGER


           On entry,  N specifies the order of the matrix C.  N must be


           at least zero.


           Unchanged on exit.

K

          K is INTEGER


           On entry with  TRANS = 'N' or 'n',  K  specifies  the number


           of  columns   of  the   matrix   A,   and  on   entry   with


           TRANS = 'C' or 'c',  K  specifies  the number of rows of the


           matrix A.  K must be at least zero.


           Unchanged on exit.

ALPHA

          ALPHA is DOUBLE PRECISION


           On entry, ALPHA specifies the scalar alpha.


           Unchanged on exit.

A

          A is COMPLEX*16 array, dimension (LDA,ka)


           where KA


           is K  when TRANS = 'N' or 'n', and is N otherwise. Before


           entry with TRANS = 'N' or 'n', the leading N--by--K part of


           the array A must contain the matrix A, otherwise the leading


           K--by--N part of the array A must contain the matrix A.


           Unchanged on exit.

LDA

          LDA is INTEGER


           On entry, LDA specifies the first dimension of A as declared


           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'


           then  LDA must be at least  max( 1, n ), otherwise  LDA must


           be at least  max( 1, k ).


           Unchanged on exit.

BETA

          BETA is DOUBLE PRECISION


           On entry, BETA specifies the scalar beta.


           Unchanged on exit.

C

          C is COMPLEX*16 array, dimension (N*(N+1)/2)


           On entry, the matrix A in RFP Format. RFP Format is


           described by TRANSR, UPLO and N. Note that the imaginary


           parts of the diagonal elements need not be set, they are


           assumed to be zero, and on exit they are set to zero.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHPCON estimates the reciprocal of the condition number of a complex


 Hermitian packed matrix A using the factorization A = U*D*U**H or


 A = L*D*L**H computed by ZHPTRF.


 An estimate is obtained for norm(inv(A)), and the reciprocal of the


 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

UPLO

          UPLO is CHARACTER*1


          Specifies whether the details of the factorization are stored


          as an upper or lower triangular matrix.


          = 'U':  Upper triangular, form is A = U*D*U**H;


          = 'L':  Lower triangular, form is A = L*D*L**H.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The block diagonal matrix D and the multipliers used to


          obtain the factor U or L as computed by ZHPTRF, stored as a


          packed triangular matrix.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZHPTRF.

ANORM

          ANORM is DOUBLE PRECISION


          The 1-norm of the original matrix A.

RCOND

          RCOND is DOUBLE PRECISION


          The reciprocal of the condition number of the matrix A,


          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an


          estimate of the 1-norm of inv(A) computed in this routine.

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHPGST reduces a complex Hermitian-definite generalized


 eigenproblem to standard form, using packed storage.


 If ITYPE = 1, the problem is A*x = lambda*B*x,


 and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)


 If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or


 B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.


 B must have been previously factorized as U**H*U or L*L**H by ZPPTRF.

ITYPE

          ITYPE is INTEGER


          = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);


          = 2 or 3: compute U*A*U**H or L**H*A*L.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored and B is factored as


                  U**H*U;


          = 'L':  Lower triangle of A is stored and B is factored as


                  L*L**H.

N

          N is INTEGER


          The order of the matrices A and B.  N >= 0.

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the upper or lower triangle of the Hermitian matrix


          A, packed columnwise in a linear array.  The j-th column of A


          is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.


          On exit, if INFO = 0, the transformed matrix, stored in the


          same format as A.

BP

          BP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The triangular factor from the Cholesky factorization of B,


          stored in the same format as A, as returned by ZPPTRF.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHPRFS improves the computed solution to a system of linear


 equations when the coefficient matrix is Hermitian indefinite


 and packed, and provides error bounds and backward error estimates


 for the solution.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

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 B and X.  NRHS >= 0.

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The upper or lower triangle of the Hermitian matrix A, packed


          columnwise in a linear array.  The j-th column of A is stored


          in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

AFP

          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The factored form of the matrix A.  AFP contains the block


          diagonal matrix D and the multipliers used to obtain the


          factor U or L from the factorization A = U*D*U**H or


          A = L*D*L**H as computed by ZHPTRF, stored as a packed


          triangular matrix.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZHPTRF.

B

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


          The right hand side matrix B.

LDB

          LDB is INTEGER


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

X

          X is COMPLEX*16 array, dimension (LDX,NRHS)


          On entry, the solution matrix X, as computed by ZHPTRS.


          On exit, the improved solution matrix X.

LDX

          LDX is INTEGER


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

FERR

          FERR is DOUBLE PRECISION array, dimension (NRHS)


          The estimated forward error bound for each solution vector


          X(j) (the j-th column of the solution matrix X).


          If XTRUE is the true solution corresponding to X(j), FERR(j)


          is an estimated upper bound for the magnitude of the largest


          element in (X(j) - XTRUE) divided by the magnitude of the


          largest element in X(j).  The estimate is as reliable as


          the estimate for RCOND, and is almost always a slight


          overestimate of the true error.

BERR

          BERR is DOUBLE PRECISION array, dimension (NRHS)


          The componentwise relative backward error of each solution


          vector X(j) (i.e., the smallest relative change in


          any element of A or B that makes X(j) an exact solution).

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

  ITMAX is the maximum number of steps of iterative refinement.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHPTRD reduces a complex Hermitian matrix A stored in packed form to


 real symmetric tridiagonal form T by a unitary similarity


 transformation: Q**H * A * Q = T.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the upper or lower triangle of the Hermitian matrix


          A, packed columnwise in a linear array.  The j-th column of A


          is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.


          On exit, if UPLO = 'U', the diagonal and first superdiagonal


          of A are overwritten by the corresponding elements of the


          tridiagonal matrix T, and the elements above the first


          superdiagonal, with the array TAU, represent the unitary


          matrix Q as a product of elementary reflectors; if UPLO


          = 'L', the diagonal and first subdiagonal of A are over-


          written by the corresponding elements of the tridiagonal


          matrix T, and the elements below the first subdiagonal, with


          the array TAU, represent the unitary matrix Q as a product


          of elementary reflectors. See Further Details.

D

          D is DOUBLE PRECISION array, dimension (N)


          The diagonal elements of the tridiagonal matrix T:


          D(i) = A(i,i).

E

          E is DOUBLE PRECISION array, dimension (N-1)


          The off-diagonal elements of the tridiagonal matrix T:


          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.

TAU

          TAU is COMPLEX*16 array, dimension (N-1)


          The scalar factors of the elementary reflectors (see Further


          Details).

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


  reflectors


     Q = H(n-1) . . . H(2) H(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+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,


  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).


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


  reflectors


     Q = H(1) H(2) . . . H(n-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(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,


  overwriting A(i+2:n,i), and tau is stored in TAU(i).

 ZHPTRF computes the factorization of a complex Hermitian packed


 matrix A using the Bunch-Kaufman diagonal pivoting method:


    A = U*D*U**H  or  A = L*D*L**H


 where U (or L) is a product of permutation and unit upper (lower)


 triangular matrices, and D is Hermitian and block diagonal with


 1-by-1 and 2-by-2 diagonal blocks.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the upper or lower triangle of the Hermitian matrix


          A, packed columnwise in a linear array.  The j-th column of A


          is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.


          On exit, the block diagonal matrix D and the multipliers used


          to obtain the factor U or L, stored as a packed triangular


          matrix overwriting A (see below for further details).

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D.


          If IPIV(k) > 0, then rows and columns k and IPIV(k) were


          interchanged and D(k,k) is a 1-by-1 diagonal block.


          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and


          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)


          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =


          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were


          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.

INFO

          INFO is INTEGER


          = 0: successful exit


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


          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization


               has been completed, but the block diagonal matrix D is


               exactly singular, and division by zero will occur if it


               is used to solve a system of equations.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  If UPLO = 'U', then A = U*D*U**H, where


     U = P(n)*U(n)* ... *P(k)U(k)* ...,


  i.e., U is a product of terms P(k)*U(k), where k decreases from n to


  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1


  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as


  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such


  that if the diagonal block D(k) is of order s (s = 1 or 2), then


             (   I    v    0   )   k-s


     U(k) =  (   0    I    0   )   s


             (   0    0    I   )   n-k


                k-s   s   n-k


  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).


  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),


  and A(k,k), and v overwrites A(1:k-2,k-1:k).


  If UPLO = 'L', then A = L*D*L**H, where


     L = P(1)*L(1)* ... *P(k)*L(k)* ...,


  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to


  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1


  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as


  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such


  that if the diagonal block D(k) is of order s (s = 1 or 2), then


             (   I    0     0   )  k-1


     L(k) =  (   0    I     0   )  s


             (   0    v     I   )  n-k-s+1


                k-1   s  n-k-s+1


  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).


  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),


  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).

J. Lewis, Boeing Computer Services Company

 ZHPTRI computes the inverse of a complex Hermitian indefinite matrix


 A in packed storage using the factorization A = U*D*U**H or


 A = L*D*L**H computed by ZHPTRF.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the details of the factorization are stored


          as an upper or lower triangular matrix.


          = 'U':  Upper triangular, form is A = U*D*U**H;


          = 'L':  Lower triangular, form is A = L*D*L**H.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the block diagonal matrix D and the multipliers


          used to obtain the factor U or L as computed by ZHPTRF,


          stored as a packed triangular matrix.


          On exit, if INFO = 0, the (Hermitian) inverse of the original


          matrix, stored as a packed triangular matrix. The j-th column


          of inv(A) is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;


          if UPLO = 'L',


             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZHPTRF.

WORK

          WORK is COMPLEX*16 array, dimension (N)

INFO

          INFO is INTEGER


          = 0: successful exit


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


          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its


               inverse could not be computed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHPTRS solves a system of linear equations A*X = B with a complex


 Hermitian matrix A stored in packed format using the factorization


 A = U*D*U**H or A = L*D*L**H computed by ZHPTRF.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the details of the factorization are stored


          as an upper or lower triangular matrix.


          = 'U':  Upper triangular, form is A = U*D*U**H;


          = 'L':  Lower triangular, form is A = L*D*L**H.

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 matrix B.  NRHS >= 0.

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The block diagonal matrix D and the multipliers used to


          obtain the factor U or L as computed by ZHPTRF, stored as a


          packed triangular matrix.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZHPTRF.

B

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


          On entry, the right hand side matrix B.


          On exit, the solution matrix X.

LDB

          LDB is INTEGER


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

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZHSEIN uses inverse iteration to find specified right and/or left


 eigenvectors of a complex upper Hessenberg matrix H.


 The right eigenvector x and the left eigenvector y of the matrix H


 corresponding to an eigenvalue w are defined by:


              H * x = w * x,     y**h * H = w * y**h


 where y**h denotes the conjugate transpose of the vector y.

SIDE

          SIDE is CHARACTER*1


          = 'R': compute right eigenvectors only;


          = 'L': compute left eigenvectors only;


          = 'B': compute both right and left eigenvectors.

EIGSRC

          EIGSRC is CHARACTER*1


          Specifies the source of eigenvalues supplied in W:


          = 'Q': the eigenvalues were found using ZHSEQR; thus, if


                 H has zero subdiagonal elements, and so is


                 block-triangular, then the j-th eigenvalue can be


                 assumed to be an eigenvalue of the block containing


                 the j-th row/column.  This property allows ZHSEIN to


                 perform inverse iteration on just one diagonal block.


          = 'N': no assumptions are made on the correspondence


                 between eigenvalues and diagonal blocks.  In this


                 case, ZHSEIN must always perform inverse iteration


                 using the whole matrix H.

INITV

          INITV is CHARACTER*1


          = 'N': no initial vectors are supplied;


          = 'U': user-supplied initial vectors are stored in the arrays


                 VL and/or VR.

SELECT

          SELECT is LOGICAL array, dimension (N)


          Specifies the eigenvectors to be computed. To select the


          eigenvector corresponding to the eigenvalue W(j),


          SELECT(j) must be set to .TRUE..

N

          N is INTEGER


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

H

          H is COMPLEX*16 array, dimension (LDH,N)


          The upper Hessenberg matrix H.


          If a NaN is detected in H, the routine will return with INFO=-6.

LDH

          LDH is INTEGER


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

W

          W is COMPLEX*16 array, dimension (N)


          On entry, the eigenvalues of H.


          On exit, the real parts of W may have been altered since


          close eigenvalues are perturbed slightly in searching for


          independent eigenvectors.

VL

          VL is COMPLEX*16 array, dimension (LDVL,MM)


          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must


          contain starting vectors for the inverse iteration for the


          left eigenvectors; the starting vector for each eigenvector


          must be in the same column in which the eigenvector will be


          stored.


          On exit, if SIDE = 'L' or 'B', the left eigenvectors


          specified by SELECT will be stored consecutively in the


          columns of VL, in the same order as their eigenvalues.


          If SIDE = 'R', VL is not referenced.

LDVL

          LDVL is INTEGER


          The leading dimension of the array VL.


          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.

VR

          VR is COMPLEX*16 array, dimension (LDVR,MM)


          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must


          contain starting vectors for the inverse iteration for the


          right eigenvectors; the starting vector for each eigenvector


          must be in the same column in which the eigenvector will be


          stored.


          On exit, if SIDE = 'R' or 'B', the right eigenvectors


          specified by SELECT will be stored consecutively in the


          columns of VR, in the same order as their eigenvalues.


          If SIDE = 'L', VR is not referenced.

LDVR

          LDVR is INTEGER


          The leading dimension of the array VR.


          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.

MM

          MM is INTEGER


          The number of columns in the arrays VL and/or VR. MM >= M.

M

          M is INTEGER


          The number of columns in the arrays VL and/or VR required to


          store the eigenvectors (= the number of .TRUE. elements in


          SELECT).

WORK

          WORK is COMPLEX*16 array, dimension (N*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

IFAILL

          IFAILL is INTEGER array, dimension (MM)


          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left


          eigenvector in the i-th column of VL (corresponding to the


          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the


          eigenvector converged satisfactorily.


          If SIDE = 'R', IFAILL is not referenced.

IFAILR

          IFAILR is INTEGER array, dimension (MM)


          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right


          eigenvector in the i-th column of VR (corresponding to the


          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the


          eigenvector converged satisfactorily.


          If SIDE = 'L', IFAILR is not referenced.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, i is the number of eigenvectors which


                failed to converge; see IFAILL and IFAILR for further


                details.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Each eigenvector is normalized so that the element of largest


  magnitude has magnitude 1; here the magnitude of a complex number


  (x,y) is taken to be |x|+|y|.

    ZHSEQR 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)*T*(QZ)**H.

JOB

          JOB is CHARACTER*1


           = 'E':  compute eigenvalues only;


           = 'S':  compute eigenvalues and the Schur form T.

COMPZ

          COMPZ is CHARACTER*1


           = 'N':  no Schur vectors are computed;


           = 'I':  Z is initialized to the unit matrix and the matrix Z


                   of Schur vectors of H is returned;


           = 'V':  Z must contain an unitary matrix Q on entry, and


                   the product Q*Z is returned.

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. ILO and IHI are normally


           set by a previous call to ZGEBAL, and then passed to ZGEHRD


           when the matrix output by ZGEBAL is reduced to Hessenberg


           form. Otherwise ILO and IHI should be set to 1 and N


           respectively.  If N > 0, then 1 <= ILO <= IHI <= N.


           If N = 0, then ILO = 1 and IHI = 0.

H

          H is COMPLEX*16 array, dimension (LDH,N)


           On entry, the upper Hessenberg matrix H.


           On exit, if INFO = 0 and JOB = 'S', H contains the upper


           triangular matrix T from the Schur decomposition (the


           Schur form). If INFO = 0 and JOB = 'E', the contents of


           H are unspecified on exit.  (The output value of H when


           INFO > 0 is given under the description of INFO below.)


           Unlike earlier versions of ZHSEQR, this subroutine may


           explicitly H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1


           or j = IHI+1, IHI+2, ... N.

LDH

          LDH is INTEGER


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

W

          W is COMPLEX*16 array, dimension (N)


           The computed eigenvalues. If JOB = 'S', 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).

Z

          Z is COMPLEX*16 array, dimension (LDZ,N)


           If COMPZ = 'N', Z is not referenced.


           If COMPZ = 'I', on entry Z need not be set and on exit,


           if INFO = 0, Z contains the unitary matrix Z of the Schur


           vectors of H.  If COMPZ = 'V', on entry Z must contain an


           N-by-N matrix Q, which is assumed to be equal to the unit


           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,


           if INFO = 0, Z contains Q*Z.


           Normally Q is the unitary matrix generated by ZUNGHR


           after the call to ZGEHRD which formed the Hessenberg matrix


           H. (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 COMPZ = 'I' or


           COMPZ = 'V', then LDZ >= MAX(1,N).  Otherwise, LDZ >= 1.

WORK

          WORK is COMPLEX*16 array, dimension (LWORK)


           On exit, if INFO = 0, 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 and delivers very good and sometimes


           optimal performance.  However, LWORK as large as 11*N


           may be required for optimal performance.  A workspace


           query is recommended to determine the optimal workspace


           size.


           If LWORK = -1, then ZHSEQR does a workspace query.


           In this case, ZHSEQR 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, the i-th argument had an illegal


                    value


             > 0:  if INFO = i, ZHSEQR failed to compute all of


                the eigenvalues.  Elements 1:ilo-1 and i+1:n of W


                contain those eigenvalues which have been


                successfully computed.  (Failures are rare.)


                If INFO > 0 and JOB = 'E', 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 JOB   = 'S', 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 COMPZ = 'V', then on exit


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


                where U is the unitary matrix in (*) (regard-


                less of the value of JOB.)


                If INFO > 0 and COMPZ = 'I', then on exit


                      (final value of Z)  = U


                where U is the unitary matrix in (*) (regard-


                less of the value of JOB.)


                If INFO > 0 and COMPZ = 'N', then Z is not


                accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

             Default values supplied by


             ILAENV(ISPEC,'ZHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).


             It is suggested that these defaults be adjusted in order


             to attain best performance in each particular


             computational environment.


            ISPEC=12: The ZLAHQR vs ZLAQR0 crossover point.


                      Default: 75. (Must be at least 11.)


            ISPEC=13: Recommended deflation window size.


                      This depends on ILO, IHI and NS.  NS is the


                      number of simultaneous shifts returned


                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)


                      The default for (IHI-ILO+1) <= 500 is NS.


                      The default for (IHI-ILO+1) >  500 is 3*NS/2.


            ISPEC=14: Nibble crossover point. (See IPARMQ for


                      details.)  Default: 14% of deflation window


                      size.


            ISPEC=15: Number of simultaneous shifts in a multishift


                      QR iteration.


                      If IHI-ILO+1 is ...


                      greater than      ...but less    ... the


                      or equal to ...      than        default is


                           1               30          NS =   2(+)


                          30               60          NS =   4(+)


                          60              150          NS =  10(+)


                         150              590          NS =  **


                         590             3000          NS =  64


                        3000             6000          NS = 128


                        6000             infinity      NS = 256


                  (+)  By default some or all matrices of this order


                       are passed to the implicit double shift routine


                       ZLAHQR and this parameter is ignored.  See


                       ISPEC=12 above and comments in IPARMQ for


                       details.


                 (**)  The asterisks (**) indicate an ad-hoc


                       function of N increasing from 10 to 64.


            ISPEC=16: Select structured matrix multiply.


                      If the number of simultaneous shifts (specified


                      by ISPEC=15) is less than 14, then the default


                      for ISPEC=16 is 0.  Otherwise the default for


                      ISPEC=16 is 2.

  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.

    ZLA_LIN_BERR computes componentwise relative backward error from


    the formula


        max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )


    where abs(Z) is the componentwise absolute value of the matrix


    or vector Z.

N

          N is INTEGER


     The number of linear equations, i.e., the order of the


     matrix A.  N >= 0.

NZ

          NZ is INTEGER


     We add (NZ+1)*SLAMCH( 'Safe minimum' ) to R(i) in the numerator to


     guard against spuriously zero residuals. Default value is N.

NRHS

          NRHS is INTEGER


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


     of the matrices AYB, RES, and BERR.  NRHS >= 0.

RES

          RES is COMPLEX*16 array, dimension (N,NRHS)


     The residual matrix, i.e., the matrix R in the relative backward


     error formula above.

AYB

          AYB is DOUBLE PRECISION array, dimension (N, NRHS)


     The denominator in the relative backward error formula above, i.e.,


     the matrix abs(op(A_s))*abs(Y) + abs(B_s). The matrices A, Y, and B


     are from iterative refinement (see zla_gerfsx_extended.f).

BERR

          BERR is DOUBLE PRECISION array, dimension (NRHS)


     The componentwise relative backward error from the formula above.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

    ZLA_WWADDW adds a vector W into a doubled-single vector (X, Y).


    This works for all extant IBM's hex and binary floating point


    arithmetic, but not for decimal.

N

          N is INTEGER


            The length of vectors X, Y, and W.

X

          X is COMPLEX*16 array, dimension (N)


            The first part of the doubled-single accumulation vector.

Y

          Y is COMPLEX*16 array, dimension (N)


            The second part of the doubled-single accumulation vector.

W

          W is COMPLEX*16 array, dimension (N)


            The vector to be added.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 Using the divide and conquer method, ZLAED0 computes all eigenvalues


 of a symmetric tridiagonal matrix which is one diagonal block of


 those from reducing a dense or band Hermitian matrix and


 corresponding eigenvectors of the dense or band matrix.

QSIZ

          QSIZ is INTEGER


         The dimension of the unitary matrix used to reduce


         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.

N

          N is INTEGER


         The dimension of the symmetric tridiagonal matrix.  N >= 0.

D

          D is DOUBLE PRECISION array, dimension (N)


         On entry, the diagonal elements of the tridiagonal matrix.


         On exit, the eigenvalues in ascending order.

E

          E is DOUBLE PRECISION array, dimension (N-1)


         On entry, the off-diagonal elements of the tridiagonal matrix.


         On exit, E has been destroyed.

Q

          Q is COMPLEX*16 array, dimension (LDQ,N)


         On entry, Q must contain an QSIZ x N matrix whose columns


         unitarily orthonormal. It is a part of the unitary matrix


         that reduces the full dense Hermitian matrix to a


         (reducible) symmetric tridiagonal matrix.

LDQ

          LDQ is INTEGER


         The leading dimension of the array Q.  LDQ >= max(1,N).

IWORK

          IWORK is INTEGER array,


         the dimension of IWORK must be at least


                      6 + 6*N + 5*N*lg N


                      ( lg( N ) = smallest integer k


                                  such that 2^k >= N )

RWORK

          RWORK is DOUBLE PRECISION array,


                               dimension (1 + 3*N + 2*N*lg N + 3*N**2)


                        ( lg( N ) = smallest integer k


                                    such that 2^k >= N )

QSTORE

          QSTORE is COMPLEX*16 array, dimension (LDQS, N)


         Used to store parts of


         the eigenvector matrix when the updating matrix multiplies


         take place.

LDQS

          LDQS is INTEGER


         The leading dimension of the array QSTORE.


         LDQS >= max(1,N).

INFO

          INFO is INTEGER


          = 0:  successful exit.


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


          > 0:  The algorithm failed to compute an eigenvalue while


                working on the submatrix lying in rows and columns


                INFO/(N+1) through mod(INFO,N+1).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZLAED7 computes the updated eigensystem of a diagonal


 matrix after modification by a rank-one symmetric matrix. This


 routine is used only for the eigenproblem which requires all


 eigenvalues and optionally eigenvectors of a dense or banded


 Hermitian matrix that has been reduced to tridiagonal form.


   T = Q(in) ( D(in) + RHO * Z*Z**H ) Q**H(in) = Q(out) * D(out) * Q**H(out)


   where Z = Q**Hu, u is a vector of length N with ones in the


   CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.


    The eigenvectors of the original matrix are stored in Q, and the


    eigenvalues are in D.  The algorithm consists of three stages:


       The first stage consists of deflating the size of the problem


       when there are multiple eigenvalues or if there is a zero in


       the Z vector.  For each such occurrence the dimension of the


       secular equation problem is reduced by one.  This stage is


       performed by the routine DLAED2.


       The second stage consists of calculating the updated


       eigenvalues. This is done by finding the roots of the secular


       equation via the routine DLAED4 (as called by SLAED3).


       This routine also calculates the eigenvectors of the current


       problem.


       The final stage consists of computing the updated eigenvectors


       directly using the updated eigenvalues.  The eigenvectors for


       the current problem are multiplied with the eigenvectors from


       the overall problem.

N

          N is INTEGER


         The dimension of the symmetric tridiagonal matrix.  N >= 0.

CUTPNT

          CUTPNT is INTEGER


         Contains the location of the last eigenvalue in the leading


         sub-matrix.  min(1,N) <= CUTPNT <= N.

QSIZ

          QSIZ is INTEGER


         The dimension of the unitary matrix used to reduce


         the full matrix to tridiagonal form.  QSIZ >= N.

TLVLS

          TLVLS is INTEGER


         The total number of merging levels in the overall divide and


         conquer tree.

CURLVL

          CURLVL is INTEGER


         The current level in the overall merge routine,


         0 <= curlvl <= tlvls.

CURPBM

          CURPBM is INTEGER


         The current problem in the current level in the overall


         merge routine (counting from upper left to lower right).

D

          D is DOUBLE PRECISION array, dimension (N)


         On entry, the eigenvalues of the rank-1-perturbed matrix.


         On exit, the eigenvalues of the repaired matrix.

Q

          Q is COMPLEX*16 array, dimension (LDQ,N)


         On entry, the eigenvectors of the rank-1-perturbed matrix.


         On exit, the eigenvectors of the repaired tridiagonal matrix.

LDQ

          LDQ is INTEGER


         The leading dimension of the array Q.  LDQ >= max(1,N).

RHO

          RHO is DOUBLE PRECISION


         Contains the subdiagonal element used to create the rank-1


         modification.

INDXQ

          INDXQ is INTEGER array, dimension (N)


         This contains the permutation which will reintegrate the


         subproblem just solved back into sorted order,


         ie. D( INDXQ( I = 1, N ) ) will be in ascending order.

IWORK

          IWORK is INTEGER array, dimension (4*N)

RWORK

          RWORK is DOUBLE PRECISION array,


                                 dimension (3*N+2*QSIZ*N)

WORK

          WORK is COMPLEX*16 array, dimension (QSIZ*N)

QSTORE

          QSTORE is DOUBLE PRECISION array, dimension (N**2+1)


         Stores eigenvectors of submatrices encountered during


         divide and conquer, packed together. QPTR points to


         beginning of the submatrices.

QPTR

          QPTR is INTEGER array, dimension (N+2)


         List of indices pointing to beginning of submatrices stored


         in QSTORE. The submatrices are numbered starting at the


         bottom left of the divide and conquer tree, from left to


         right and bottom to top.

PRMPTR

          PRMPTR is INTEGER array, dimension (N lg N)


         Contains a list of pointers which indicate where in PERM a


         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)


         indicates the size of the permutation and also the size of


         the full, non-deflated problem.

PERM

          PERM is INTEGER array, dimension (N lg N)


         Contains the permutations (from deflation and sorting) to be


         applied to each eigenblock.

GIVPTR

          GIVPTR is INTEGER array, dimension (N lg N)


         Contains a list of pointers which indicate where in GIVCOL a


         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)


         indicates the number of Givens rotations.

GIVCOL

          GIVCOL is INTEGER array, dimension (2, N lg N)


         Each pair of numbers indicates a pair of columns to take place


         in a Givens rotation.

GIVNUM

          GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)


         Each number indicates the S value to be used in the


         corresponding Givens rotation.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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


          > 0:  if INFO = 1, an eigenvalue did not converge

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZLAED8 merges the two sets of eigenvalues together into a single


 sorted set.  Then it tries to deflate the size of the problem.


 There are two ways in which deflation can occur:  when two or more


 eigenvalues are close together or if there is a tiny element in the


 Z vector.  For each such occurrence the order of the related secular


 equation problem is reduced by one.

K

          K is INTEGER


         Contains the number of non-deflated eigenvalues.


         This is the order of the related secular equation.

N

          N is INTEGER


         The dimension of the symmetric tridiagonal matrix.  N >= 0.

QSIZ

          QSIZ is INTEGER


         The dimension of the unitary matrix used to reduce


         the dense or band matrix to tridiagonal form.


         QSIZ >= N if ICOMPQ = 1.

Q

          Q is COMPLEX*16 array, dimension (LDQ,N)


         On entry, Q contains the eigenvectors of the partially solved


         system which has been previously updated in matrix


         multiplies with other partially solved eigensystems.


         On exit, Q contains the trailing (N-K) updated eigenvectors


         (those which were deflated) in its last N-K columns.

LDQ

          LDQ is INTEGER


         The leading dimension of the array Q.  LDQ >= max( 1, N ).

D

          D is DOUBLE PRECISION array, dimension (N)


         On entry, D contains the eigenvalues of the two submatrices to


         be combined.  On exit, D contains the trailing (N-K) updated


         eigenvalues (those which were deflated) sorted into increasing


         order.

RHO

          RHO is DOUBLE PRECISION


         Contains the off diagonal element associated with the rank-1


         cut which originally split the two submatrices which are now


         being recombined. RHO is modified during the computation to


         the value required by DLAED3.

CUTPNT

          CUTPNT is INTEGER


         Contains the location of the last eigenvalue in the leading


         sub-matrix.  MIN(1,N) <= CUTPNT <= N.

Z

          Z is DOUBLE PRECISION array, dimension (N)


         On input this vector contains the updating vector (the last


         row of the first sub-eigenvector matrix and the first row of


         the second sub-eigenvector matrix).  The contents of Z are


         destroyed during the updating process.

DLAMDA

          DLAMDA is DOUBLE PRECISION array, dimension (N)


         Contains a copy of the first K eigenvalues which will be used


         by DLAED3 to form the secular equation.

Q2

          Q2 is COMPLEX*16 array, dimension (LDQ2,N)


         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,


         Contains a copy of the first K eigenvectors which will be used


         by DLAED7 in a matrix multiply (DGEMM) to update the new


         eigenvectors.

LDQ2

          LDQ2 is INTEGER


         The leading dimension of the array Q2.  LDQ2 >= max( 1, N ).

W

          W is DOUBLE PRECISION array, dimension (N)


         This will hold the first k values of the final


         deflation-altered z-vector and will be passed to DLAED3.

INDXP

          INDXP is INTEGER array, dimension (N)


         This will contain the permutation used to place deflated


         values of D at the end of the array. On output INDXP(1:K)


         points to the nondeflated D-values and INDXP(K+1:N)


         points to the deflated eigenvalues.

INDX

          INDX is INTEGER array, dimension (N)


         This will contain the permutation used to sort the contents of


         D into ascending order.

INDXQ

          INDXQ is INTEGER array, dimension (N)


         This contains the permutation which separately sorts the two


         sub-problems in D into ascending order.  Note that elements in


         the second half of this permutation must first have CUTPNT


         added to their values in order to be accurate.

PERM

          PERM is INTEGER array, dimension (N)


         Contains the permutations (from deflation and sorting) to be


         applied to each eigenblock.

GIVPTR

          GIVPTR is INTEGER


         Contains the number of Givens rotations which took place in


         this subproblem.

GIVCOL

          GIVCOL is INTEGER array, dimension (2, N)


         Each pair of numbers indicates a pair of columns to take place


         in a Givens rotation.

GIVNUM

          GIVNUM is DOUBLE PRECISION array, dimension (2, N)


         Each number indicates the S value to be used in the


         corresponding Givens rotation.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZLALS0 applies back the multiplying factors of either the left or the


 right singular vector matrix of a diagonal matrix appended by a row


 to the right hand side matrix B in solving the least squares problem


 using the divide-and-conquer SVD approach.


 For the left singular vector matrix, three types of orthogonal


 matrices are involved:


 (1L) Givens rotations: the number of such rotations is GIVPTR; the


      pairs of columns/rows they were applied to are stored in GIVCOL;


      and the C- and S-values of these rotations are stored in GIVNUM.


 (2L) Permutation. The (NL+1)-st row of B is to be moved to the first


      row, and for J=2:N, PERM(J)-th row of B is to be moved to the


      J-th row.


 (3L) The left singular vector matrix of the remaining matrix.


 For the right singular vector matrix, four types of orthogonal


 matrices are involved:


 (1R) The right singular vector matrix of the remaining matrix.


 (2R) If SQRE = 1, one extra Givens rotation to generate the right


      null space.


 (3R) The inverse transformation of (2L).


 (4R) The inverse transformation of (1L).

ICOMPQ

          ICOMPQ is INTEGER


         Specifies whether singular vectors are to be computed in


         factored form:


         = 0: Left singular vector matrix.


         = 1: Right singular vector matrix.

NL

          NL is INTEGER


         The row dimension of the upper block. NL >= 1.

NR

          NR is INTEGER


         The row dimension of the lower block. NR >= 1.

SQRE

          SQRE is INTEGER


         = 0: the lower block is an NR-by-NR square matrix.


         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.


         The bidiagonal matrix has row dimension N = NL + NR + 1,


         and column dimension M = N + SQRE.

NRHS

          NRHS is INTEGER


         The number of columns of B and BX. NRHS must be at least 1.

B

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


         On input, B contains the right hand sides of the least


         squares problem in rows 1 through M. On output, B contains


         the solution X in rows 1 through N.

LDB

          LDB is INTEGER


         The leading dimension of B. LDB must be at least


         max(1,MAX( M, N ) ).

BX

          BX is COMPLEX*16 array, dimension ( LDBX, NRHS )

LDBX

          LDBX is INTEGER


         The leading dimension of BX.

PERM

          PERM is INTEGER array, dimension ( N )


         The permutations (from deflation and sorting) applied


         to the two blocks.

GIVPTR

          GIVPTR is INTEGER


         The number of Givens rotations which took place in this


         subproblem.

GIVCOL

          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )


         Each pair of numbers indicates a pair of rows/columns


         involved in a Givens rotation.

LDGCOL

          LDGCOL is INTEGER


         The leading dimension of GIVCOL, must be at least N.

GIVNUM

          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )


         Each number indicates the C or S value used in the


         corresponding Givens rotation.

LDGNUM

          LDGNUM is INTEGER


         The leading dimension of arrays DIFR, POLES and


         GIVNUM, must be at least K.

POLES

          POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )


         On entry, POLES(1:K, 1) contains the new singular


         values obtained from solving the secular equation, and


         POLES(1:K, 2) is an array containing the poles in the secular


         equation.

DIFL

          DIFL is DOUBLE PRECISION array, dimension ( K ).


         On entry, DIFL(I) is the distance between I-th updated


         (undeflated) singular value and the I-th (undeflated) old


         singular value.

DIFR

          DIFR is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).


         On entry, DIFR(I, 1) contains the distances between I-th


         updated (undeflated) singular value and the I+1-th


         (undeflated) old singular value. And DIFR(I, 2) is the


         normalizing factor for the I-th right singular vector.

Z

          Z is DOUBLE PRECISION array, dimension ( K )


         Contain the components of the deflation-adjusted updating row


         vector.

K

          K is INTEGER


         Contains the dimension of the non-deflated matrix,


         This is the order of the related secular equation. 1 <= K <=N.

C

          C is DOUBLE PRECISION


         C contains garbage if SQRE =0 and the C-value of a Givens


         rotation related to the right null space if SQRE = 1.

S

          S is DOUBLE PRECISION


         S contains garbage if SQRE =0 and the S-value of a Givens


         rotation related to the right null space if SQRE = 1.

RWORK

          RWORK is DOUBLE PRECISION array, dimension


         ( K*(1+NRHS) + 2*NRHS )

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

 ZLALSA is an itermediate step in solving the least squares problem


 by computing the SVD of the coefficient matrix in compact form (The


 singular vectors are computed as products of simple orthorgonal


 matrices.).


 If ICOMPQ = 0, ZLALSA applies the inverse of the left singular vector


 matrix of an upper bidiagonal matrix to the right hand side; and if


 ICOMPQ = 1, ZLALSA applies the right singular vector matrix to the


 right hand side. The singular vector matrices were generated in


 compact form by ZLALSA.

ICOMPQ

          ICOMPQ is INTEGER


         Specifies whether the left or the right singular vector


         matrix is involved.


         = 0: Left singular vector matrix


         = 1: Right singular vector matrix

SMLSIZ

          SMLSIZ is INTEGER


         The maximum size of the subproblems at the bottom of the


         computation tree.

N

          N is INTEGER


         The row and column dimensions of the upper bidiagonal matrix.

NRHS

          NRHS is INTEGER


         The number of columns of B and BX. NRHS must be at least 1.

B

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


         On input, B contains the right hand sides of the least


         squares problem in rows 1 through M.


         On output, B contains the solution X in rows 1 through N.

LDB

          LDB is INTEGER


         The leading dimension of B in the calling subprogram.


         LDB must be at least max(1,MAX( M, N ) ).

BX

          BX is COMPLEX*16 array, dimension ( LDBX, NRHS )


         On exit, the result of applying the left or right singular


         vector matrix to B.

LDBX

          LDBX is INTEGER


         The leading dimension of BX.

U

          U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).


         On entry, U contains the left singular vector matrices of all


         subproblems at the bottom level.

LDU

          LDU is INTEGER, LDU = > N.


         The leading dimension of arrays U, VT, DIFL, DIFR,


         POLES, GIVNUM, and Z.

VT

          VT is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).


         On entry, VT**H contains the right singular vector matrices of


         all subproblems at the bottom level.

K

          K is INTEGER array, dimension ( N ).

DIFL

          DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ).


         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.

DIFR

          DIFR is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).


         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record


         distances between singular values on the I-th level and


         singular values on the (I -1)-th level, and DIFR(*, 2 * I)


         record the normalizing factors of the right singular vectors


         matrices of subproblems on I-th level.

Z

          Z is DOUBLE PRECISION array, dimension ( LDU, NLVL ).


         On entry, Z(1, I) contains the components of the deflation-


         adjusted updating row vector for subproblems on the I-th


         level.

POLES

          POLES is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).


         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old


         singular values involved in the secular equations on the I-th


         level.

GIVPTR

          GIVPTR is INTEGER array, dimension ( N ).


         On entry, GIVPTR( I ) records the number of Givens


         rotations performed on the I-th problem on the computation


         tree.

GIVCOL

          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).


         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the


         locations of Givens rotations performed on the I-th level on


         the computation tree.

LDGCOL

          LDGCOL is INTEGER, LDGCOL = > N.


         The leading dimension of arrays GIVCOL and PERM.

PERM

          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).


         On entry, PERM(*, I) records permutations done on the I-th


         level of the computation tree.

GIVNUM

          GIVNUM is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).


         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-


         values of Givens rotations performed on the I-th level on the


         computation tree.

C

          C is DOUBLE PRECISION array, dimension ( N ).


         On entry, if the I-th subproblem is not square,


         C( I ) contains the C-value of a Givens rotation related to


         the right null space of the I-th subproblem.

S

          S is DOUBLE PRECISION array, dimension ( N ).


         On entry, if the I-th subproblem is not square,


         S( I ) contains the S-value of a Givens rotation related to


         the right null space of the I-th subproblem.

RWORK

          RWORK is DOUBLE PRECISION array, dimension at least


         MAX( (SMLSZ+1)*NRHS*3, N*(1+NRHS) + 2*NRHS ).

IWORK

          IWORK is INTEGER array, dimension (3*N)

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

 ZLALSD uses the singular value decomposition of A to solve the least


 squares problem of finding X to minimize the Euclidean norm of each


 column of A*X-B, where A is N-by-N upper bidiagonal, and X and B


 are N-by-NRHS. The solution X overwrites B.


 The singular values of A smaller than RCOND times the largest


 singular value are treated as zero in solving the least squares


 problem; in this case a minimum norm solution is returned.


 The actual singular values are returned in D in ascending order.


 This code makes very mild assumptions about floating point


 arithmetic. It will work on machines with a guard digit in


 add/subtract, or on those binary machines without guard digits


 which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.


 It could conceivably fail on hexadecimal or decimal machines


 without guard digits, but we know of none.

UPLO

          UPLO is CHARACTER*1


         = 'U': D and E define an upper bidiagonal matrix.


         = 'L': D and E define a  lower bidiagonal matrix.

SMLSIZ

          SMLSIZ is INTEGER


         The maximum size of the subproblems at the bottom of the


         computation tree.

N

          N is INTEGER


         The dimension of the  bidiagonal matrix.  N >= 0.

NRHS

          NRHS is INTEGER


         The number of columns of B. NRHS must be at least 1.

D

          D is DOUBLE PRECISION array, dimension (N)


         On entry D contains the main diagonal of the bidiagonal


         matrix. On exit, if INFO = 0, D contains its singular values.

E

          E is DOUBLE PRECISION array, dimension (N-1)


         Contains the super-diagonal entries of the bidiagonal matrix.


         On exit, E has been destroyed.

B

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


         On input, B contains the right hand sides of the least


         squares problem. On output, B contains the solution X.

LDB

          LDB is INTEGER


         The leading dimension of B in the calling subprogram.


         LDB must be at least max(1,N).

RCOND

          RCOND is DOUBLE PRECISION


         The singular values of A less than or equal to RCOND times


         the largest singular value are treated as zero in solving


         the least squares problem. If RCOND is negative,


         machine precision is used instead.


         For example, if diag(S)*X=B were the least squares problem,


         where diag(S) is a diagonal matrix of singular values, the


         solution would be X(i) = B(i) / S(i) if S(i) is greater than


         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to


         RCOND*max(S).

RANK

          RANK is INTEGER


         The number of singular values of A greater than RCOND times


         the largest singular value.

WORK

          WORK is COMPLEX*16 array, dimension (N * NRHS)

RWORK

          RWORK is DOUBLE PRECISION array, dimension at least


         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +


         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),


         where


         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )

IWORK

          IWORK is INTEGER array, dimension at least


         (3*N*NLVL + 11*N).

INFO

          INFO is INTEGER


         = 0:  successful exit.


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


         > 0:  The algorithm failed to compute a singular value while


               working on the submatrix lying in rows and columns


               INFO/(N+1) through MOD(INFO,N+1).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

 ZLANHF  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 in RFP format.

ZLANHF

    ZLANHF = ( 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  matrix norm.

NORM

          NORM is CHARACTER


            Specifies the value to be returned in ZLANHF as described


            above.

TRANSR

          TRANSR is CHARACTER


            Specifies whether the RFP format of A is normal or


            conjugate-transposed format.


            = 'N':  RFP format is Normal


            = 'C':  RFP format is Conjugate-transposed

UPLO

          UPLO is CHARACTER


            On entry, UPLO specifies whether the RFP matrix A came from


            an upper or lower triangular matrix as follows:


            UPLO = 'U' or 'u' RFP A came from an upper triangular


            matrix


            UPLO = 'L' or 'l' RFP A came from a  lower triangular


            matrix

N

          N is INTEGER


            The order of the matrix A.  N >= 0.  When N = 0, ZLANHF is


            set to zero.

A

          A is COMPLEX*16 array, dimension ( N*(N+1)/2 );


            On entry, the matrix A in RFP Format.


            RFP Format is described by TRANSR, UPLO and N as follows:


            If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;


            K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If


            TRANSR = 'C' then RFP is the Conjugate-transpose of RFP A


            as defined when TRANSR = 'N'. The contents of RFP A are


            defined by UPLO as follows: If UPLO = 'U' the RFP A


            contains the ( N*(N+1)/2 ) elements of upper packed A


            either in normal or conjugate-transpose Format. If


            UPLO = 'L' the RFP A contains the ( N*(N+1) /2 ) elements


            of lower packed A either in normal or conjugate-transpose


            Format. The LDA of RFP A is (N+1)/2 when TRANSR = 'C'. When


            TRANSR is 'N' the LDA is N+1 when N is even and is N when


            is odd. See the Note below for more details.


            Unchanged on exit.

WORK

          WORK is DOUBLE PRECISION array, dimension (LWORK),


            where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,


            WORK is not referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  We first consider Standard Packed Format when N is even.


  We give an example where N = 6.


      AP is Upper             AP is Lower


   00 01 02 03 04 05       00


      11 12 13 14 15       10 11


         22 23 24 25       20 21 22


            33 34 35       30 31 32 33


               44 45       40 41 42 43 44


                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last


  three columns of AP upper. The lower triangle A(4:6,0:2) consists of


  conjugate-transpose of the first three columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:2,0:2) consists of


  conjugate-transpose of the last three columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N even and TRANSR = 'N'.


         RFP A                   RFP A


                                -- -- --


        03 04 05                33 43 53


                                   -- --


        13 14 15                00 44 54


                                      --


        23 24 25                10 11 55


        33 34 35                20 21 22


        --


        00 44 45                30 31 32


        -- --


        01 11 55                40 41 42


        -- -- --


        02 12 22                50 51 52


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- -- --                -- -- -- -- -- --


     03 13 23 33 00 01 02    33 00 10 20 30 40 50


     -- -- -- -- --                -- -- -- -- --


     04 14 24 34 44 11 12    43 44 11 21 31 41 51


     -- -- -- -- -- --                -- -- -- --


     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We next  consider Standard Packed Format when N is odd.


  We give an example where N = 5.


     AP is Upper                 AP is Lower


   00 01 02 03 04              00


      11 12 13 14              10 11


         22 23 24              20 21 22


            33 34              30 31 32 33


               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last


  three columns of AP upper. The lower triangle A(3:4,0:1) consists of


  conjugate-transpose of the first two   columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:1,1:2) consists of


  conjugate-transpose of the last two   columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N odd  and TRANSR = 'N'.


         RFP A                   RFP A


                                   -- --


        02 03 04                00 33 43


                                      --


        12 13 14                10 11 44


        22 23 24                20 21 22


        --


        00 33 34                30 31 32


        -- --


        01 11 44                40 41 42


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- --                   -- -- -- -- -- --


     02 12 22 00 01             00 10 20 30 40 50


     -- -- -- --                   -- -- -- -- --


     03 13 23 33 11             33 11 21 31 41 51


     -- -- -- -- --                   -- -- -- --


     04 14 24 34 44             43 44 22 32 42 52

 ZLARSCL2 performs a reciprocal diagonal scaling on a matrix:


   x <-- inv(D) * x


 where the DOUBLE PRECISION diagonal matrix D is stored as a vector.


 Eventually to be replaced by BLAS_zge_diag_scale in the new BLAS


 standard.

M

          M is INTEGER


     The number of rows of D and X. M >= 0.

N

          N is INTEGER


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

D

          D is DOUBLE PRECISION array, length M


     Diagonal matrix D, stored as a vector of length M.

X

          X is COMPLEX*16 array, dimension (LDX,N)


     On entry, the matrix X to be scaled by D.


     On exit, the scaled matrix.

LDX

          LDX is INTEGER


     The leading dimension of the matrix X. LDX >= M.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


 H is a product of k elementary reflectors as returned by ZTZRZF.

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.

L

          L is INTEGER


          The number of entries of the vector V containing


          the meaningful part of the Householder vectors.


          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

V

          V is COMPLEX*16 array, dimension (1+(L-1)*abs(INCV))


          The vector v in the representation of H as returned by


          ZTZRZF. V is not used if TAU = 0.

INCV

          INCV is INTEGER


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

TAU

          TAU is COMPLEX*16


          The value tau in the representation of H.

C

          C is COMPLEX*16 array, dimension (LDC,N)


          On entry, the M-by-N matrix C.


          On exit, C is overwritten by the matrix H * C if SIDE = 'L',


          or C * H if SIDE = 'R'.

LDC

          LDC is INTEGER


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

WORK

          WORK is COMPLEX*16 array, dimension


                         (N) if SIDE = 'L'


                      or (M) if SIDE = 'R'

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

 ZLARZB applies a complex block reflector H or its transpose H**H


 to a complex distributed M-by-N  C from the left or the right.


 Currently, only STOREV = 'R' and DIRECT = 'B' are supported.

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, not supported yet)


          = '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                        (not supported yet)


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

L

          L is INTEGER


          The number of columns of the matrix V containing the


          meaningful part of the Householder reflectors.


          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

V

          V is COMPLEX*16 array, dimension (LDV,NV).


          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.

T

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


          The triangular K-by-K matrix T in the representation of the


          block reflector.

LDT

          LDT is INTEGER


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

C

          C is COMPLEX*16 array, dimension (LDC,N)


          On entry, the M-by-N matrix C.


          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.

LDC

          LDC is INTEGER


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

WORK

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

LDWORK

          LDWORK is INTEGER


          The leading dimension of the array WORK.


          If SIDE = 'L', LDWORK >= max(1,N);


          if SIDE = 'R', LDWORK >= max(1,M).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

 ZLARZT 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


 Currently, only STOREV = 'R' and DIRECT = 'B' are supported.

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, not supported yet)


          = '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                        (not supported yet)


          = 'R': rowwise

N

          N is INTEGER


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

K

          K is INTEGER


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


          elementary reflectors). K >= 1.

V

          V is COMPLEX*16 array, dimension


                               (LDV,K) if STOREV = 'C'


                               (LDV,N) if STOREV = 'R'


          The matrix V. See further details.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.

TAU

          TAU is COMPLEX*16 array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i).

T

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


          The k by k triangular factor T of the block reflector.


          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is


          lower triangular. The rest of the array is not used.

LDT

          LDT is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

  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_____


         ( v1 v2 v3 )                        /            \


         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )


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


         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )


         ( v1 v2 v3 )


            .  .  .


            .  .  .


            1  .  .


               1  .


                  1


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


                                                        ______V_____


            1                                          /            \


            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )


            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )


            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )


            .  .  .


         ( v1 v2 v3 )


         ( v1 v2 v3 )


     V = ( v1 v2 v3 )


         ( v1 v2 v3 )


         ( v1 v2 v3 )

 ZLASCL2 performs a diagonal scaling on a matrix:


   x <-- D * x


 where the DOUBLE PRECISION diagonal matrix D is stored as a vector.


 Eventually to be replaced by BLAS_zge_diag_scale in the new BLAS


 standard.

M

          M is INTEGER


     The number of rows of D and X. M >= 0.

N

          N is INTEGER


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

D

          D is DOUBLE PRECISION array, length M


     Diagonal matrix D, stored as a vector of length M.

X

          X is COMPLEX*16 array, dimension (LDX,N)


     On entry, the matrix X to be scaled by D.


     On exit, the scaled matrix.

LDX

          LDX is INTEGER


     The leading dimension of the matrix X. LDX >= M.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZLATRZ factors the M-by-(M+L) complex upper trapezoidal matrix


 [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R  0 ) * Z by means


 of unitary transformations, where  Z is an (M+L)-by-(M+L) unitary


 matrix and, R and A1 are M-by-M upper triangular matrices.

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.

L

          L is INTEGER


          The number of columns of the matrix A containing the


          meaningful part of the Householder vectors. N-M >= L >= 0.

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, the leading M-by-N upper trapezoidal part of the


          array A must contain the matrix to be factorized.


          On exit, the leading M-by-M upper triangular part of A


          contains the upper triangular matrix R, and elements N-L+1 to


          N of the first M rows of A, with the array TAU, represent the


          unitary matrix Z as a product of M elementary reflectors.

LDA

          LDA is INTEGER


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

TAU

          TAU is COMPLEX*16 array, dimension (M)


          The scalar factors of the elementary reflectors.

WORK

          WORK is COMPLEX*16 array, dimension (M)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

  The factorization is obtained by Householder's method.  The kth


  transformation matrix, Z( k ), which is used to introduce zeros into


  the ( m - k + 1 )th row of A, is given in the form


     Z( k ) = ( I     0   ),


              ( 0  T( k ) )


  where


     T( k ) = I - tau*u( k )*u( k )**H,   u( k ) = (   1    ),


                                                 (   0    )


                                                 ( z( k ) )


  tau is a scalar and z( k ) is an l element vector. tau and z( k )


  are chosen to annihilate the elements of the kth row of A2.


  The scalar tau is returned in the kth element of TAU and the vector


  u( k ) in the kth row of A2, such that the elements of z( k ) are


  in  a( k, l + 1 ), ..., a( k, n ). The elements of R are returned in


  the upper triangular part of A1.


  Z is given by


     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).

 ZPBCON estimates the reciprocal of the condition number (in the


 1-norm) of a complex Hermitian positive definite band matrix using


 the Cholesky factorization A = U**H*U or A = L*L**H computed by


 ZPBTRF.


 An estimate is obtained for norm(inv(A)), and the reciprocal of the


 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangular factor stored in AB;


          = 'L':  Lower triangular factor stored in AB.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


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

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          The triangular factor U or L from the Cholesky factorization


          A = U**H*U or A = L*L**H of the band matrix A, stored in the


          first KD+1 rows of the array.  The j-th column of U or L is


          stored in the j-th column of the array AB as follows:


          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;


          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).

LDAB

          LDAB is INTEGER


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

ANORM

          ANORM is DOUBLE PRECISION


          The 1-norm (or infinity-norm) of the Hermitian band matrix A.

RCOND

          RCOND is DOUBLE PRECISION


          The reciprocal of the condition number of the matrix A,


          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an


          estimate of the 1-norm of inv(A) computed in this routine.

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPBEQU computes row and column scalings intended to equilibrate a


 Hermitian positive definite band matrix A and reduce its condition


 number (with respect to the two-norm).  S contains the scale factors,


 S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with


 elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This


 choice of S puts the condition number of B within a factor N of the


 smallest possible condition number over all possible diagonal


 scalings.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangular of A is stored;


          = 'L':  Lower triangular of A is stored.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          The upper or lower triangle of the Hermitian band matrix A,


          stored in the first 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 A.  LDAB >= KD+1.

S

          S is DOUBLE PRECISION array, dimension (N)


          If INFO = 0, S contains the scale factors for A.

SCOND

          SCOND is DOUBLE PRECISION


          If INFO = 0, S contains the ratio of the smallest S(i) to


          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too


          large nor too small, it is not worth scaling by S.

AMAX

          AMAX is DOUBLE PRECISION


          Absolute value of largest matrix element.  If AMAX is very


          close to overflow or very close to underflow, the matrix


          should be scaled.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the i-th diagonal element is nonpositive.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPBRFS improves the computed solution to a system of linear


 equations when the coefficient matrix is Hermitian positive definite


 and banded, and provides error bounds and backward error estimates


 for the solution.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

NRHS

          NRHS is INTEGER


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


          of the matrices B and X.  NRHS >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          The upper or lower triangle of the Hermitian band matrix A,


          stored in the first 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.

AFB

          AFB is COMPLEX*16 array, dimension (LDAFB,N)


          The triangular factor U or L from the Cholesky factorization


          A = U**H*U or A = L*L**H of the band matrix A as computed by


          ZPBTRF, in the same storage format as A (see AB).

LDAFB

          LDAFB is INTEGER


          The leading dimension of the array AFB.  LDAFB >= KD+1.

B

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


          The right hand side matrix B.

LDB

          LDB is INTEGER


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

X

          X is COMPLEX*16 array, dimension (LDX,NRHS)


          On entry, the solution matrix X, as computed by ZPBTRS.


          On exit, the improved solution matrix X.

LDX

          LDX is INTEGER


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

FERR

          FERR is DOUBLE PRECISION array, dimension (NRHS)


          The estimated forward error bound for each solution vector


          X(j) (the j-th column of the solution matrix X).


          If XTRUE is the true solution corresponding to X(j), FERR(j)


          is an estimated upper bound for the magnitude of the largest


          element in (X(j) - XTRUE) divided by the magnitude of the


          largest element in X(j).  The estimate is as reliable as


          the estimate for RCOND, and is almost always a slight


          overestimate of the true error.

BERR

          BERR is DOUBLE PRECISION array, dimension (NRHS)


          The componentwise relative backward error of each solution


          vector X(j) (i.e., the smallest relative change in


          any element of A or B that makes X(j) an exact solution).

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

  ITMAX is the maximum number of steps of iterative refinement.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPBSTF computes a split Cholesky factorization of a complex


 Hermitian positive definite band matrix A.


 This routine is designed to be used in conjunction with ZHBGST.


 The factorization has the form  A = S**H*S  where S is a band matrix


 of the same bandwidth as A and the following structure:


   S = ( U    )


       ( M  L )


 where U is upper triangular of order m = (n+kd)/2, and L is lower


 triangular of order n-m.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          On entry, the upper or lower triangle of the Hermitian 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 factor S from the split Cholesky


          factorization A = S**H*S. See Further Details.

LDAB

          LDAB is INTEGER


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

INFO

          INFO is INTEGER


          = 0: successful exit


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


          > 0: if INFO = i, the factorization could not be completed,


               because the updated element a(i,i) was negative; the


               matrix A is not positive definite.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The band storage scheme is illustrated by the following example, when


  N = 7, KD = 2:


  S = ( s11  s12  s13                     )


      (      s22  s23  s24                )


      (           s33  s34                )


      (                s44                )


      (           s53  s54  s55           )


      (                s64  s65  s66      )


      (                     s75  s76  s77 )


  If UPLO = 'U', the array AB holds:


  on entry:                          on exit:


   *    *   a13  a24  a35  a46  a57   *    *   s13  s24  s53**H s64**H s75**H


   *   a12  a23  a34  a45  a56  a67   *   s12  s23  s34  s54**H s65**H s76**H


  a11  a22  a33  a44  a55  a66  a77  s11  s22  s33  s44  s55    s66    s77


  If UPLO = 'L', the array AB holds:


  on entry:                          on exit:


  a11  a22  a33  a44  a55  a66  a77  s11    s22    s33    s44  s55  s66  s77


  a21  a32  a43  a54  a65  a76   *   s12**H s23**H s34**H s54  s65  s76   *


  a31  a42  a53  a64  a64   *    *   s13**H s24**H s53    s64  s75   *    *


  Array elements marked * are not used by the routine; s12**H denotes


  conjg(s12); the diagonal elements of S are real.

 ZPBTF2 computes the Cholesky factorization of a complex Hermitian


 positive definite band matrix A.


 The factorization has the form


    A = U**H * U ,  if UPLO = 'U', or


    A = L  * L**H,  if UPLO = 'L',


 where U is an upper triangular matrix, U**H is the conjugate transpose


 of U, and L is lower triangular.


 This is the unblocked version of the algorithm, calling Level 2 BLAS.

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.

KD

          KD is INTEGER


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


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

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          On entry, the upper or lower triangle of the Hermitian 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.

INFO

          INFO is INTEGER


          = 0: successful exit


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


          > 0: if INFO = k, the leading minor of order k is not


               positive definite, and the factorization could not be


               completed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The band storage scheme is illustrated by the following example, when


  N = 6, KD = 2, and UPLO = 'U':


  On entry:                       On exit:


      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46


      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56


     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66


  Similarly, if UPLO = 'L' the format of A is as follows:


  On entry:                       On exit:


     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66


     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *


     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *


  Array elements marked * are not used by the routine.

 ZPBTRF computes the Cholesky factorization of a complex Hermitian


 positive definite band matrix A.


 The factorization has the form


    A = U**H * U,  if UPLO = 'U', or


    A = L  * L**H,  if UPLO = 'L',


 where U is an upper triangular matrix and L is lower triangular.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          On entry, the upper or lower triangle of the Hermitian 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.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the leading minor of order i is not


                positive definite, and the factorization could not be


                completed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The band storage scheme is illustrated by the following example, when


  N = 6, KD = 2, and UPLO = 'U':


  On entry:                       On exit:


      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46


      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56


     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66


  Similarly, if UPLO = 'L' the format of A is as follows:


  On entry:                       On exit:


     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66


     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *


     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *


  Array elements marked * are not used by the routine.

Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989

 ZPBTRS solves a system of linear equations A*X = B with a Hermitian


 positive definite band matrix A using the Cholesky factorization


 A = U**H *U or A = L*L**H computed by ZPBTRF.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangular factor stored in AB;


          = 'L':  Lower triangular factor stored in AB.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of superdiagonals of the matrix A if UPLO = 'U',


          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

NRHS

          NRHS is INTEGER


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


          of the matrix B.  NRHS >= 0.

AB

          AB is COMPLEX*16 array, dimension (LDAB,N)


          The triangular factor U or L from the Cholesky factorization


          A = U**H *U or A = L*L**H of the band matrix A, stored in the


          first KD+1 rows of the array.  The j-th column of U or L is


          stored in the j-th column of the array AB as follows:


          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;


          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).

LDAB

          LDAB is INTEGER


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

B

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


          On entry, the right hand side matrix B.


          On exit, the solution matrix X.

LDB

          LDB is INTEGER


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

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPFTRF computes the Cholesky factorization of a complex Hermitian


 positive definite matrix A.


 The factorization has the form


    A = U**H * U,  if UPLO = 'U', or


    A = L  * L**H,  if UPLO = 'L',


 where U is an upper triangular matrix and L is lower triangular.


 This is the block version of the algorithm, calling Level 3 BLAS.

TRANSR

          TRANSR is CHARACTER*1


          = 'N':  The Normal TRANSR of RFP A is stored;


          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of RFP A is stored;


          = 'L':  Lower triangle of RFP A is stored.

N

          N is INTEGER


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

A

          A is COMPLEX*16 array, dimension ( N*(N+1)/2 );


          On entry, the Hermitian matrix A in RFP format. RFP format is


          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'


          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is


          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is


          the Conjugate-transpose of RFP A as defined when


          TRANSR = 'N'. The contents of RFP A are defined by UPLO as


          follows: If UPLO = 'U' the RFP A contains the nt elements of


          upper packed A. If UPLO = 'L' the RFP A contains the elements


          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =


          'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N


          is odd. See the Note below for more details.


          On exit, if INFO = 0, the factor U or L from the Cholesky


          factorization RFP A = U**H*U or RFP A = L*L**H.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the leading minor of order i is not


                positive definite, and the factorization could not be


                completed.


  Further Notes on RFP Format:


  ============================


  We first consider Standard Packed Format when N is even.


  We give an example where N = 6.


     AP is Upper             AP is Lower


   00 01 02 03 04 05       00


      11 12 13 14 15       10 11


         22 23 24 25       20 21 22


            33 34 35       30 31 32 33


               44 45       40 41 42 43 44


                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last


  three columns of AP upper. The lower triangle A(4:6,0:2) consists of


  conjugate-transpose of the first three columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:2,0:2) consists of


  conjugate-transpose of the last three columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N even and TRANSR = 'N'.


         RFP A                   RFP A


                                -- -- --


        03 04 05                33 43 53


                                   -- --


        13 14 15                00 44 54


                                      --


        23 24 25                10 11 55


        33 34 35                20 21 22


        --


        00 44 45                30 31 32


        -- --


        01 11 55                40 41 42


        -- -- --


        02 12 22                50 51 52


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- -- --                -- -- -- -- -- --


     03 13 23 33 00 01 02    33 00 10 20 30 40 50


     -- -- -- -- --                -- -- -- -- --


     04 14 24 34 44 11 12    43 44 11 21 31 41 51


     -- -- -- -- -- --                -- -- -- --


     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We next  consider Standard Packed Format when N is odd.


  We give an example where N = 5.


     AP is Upper                 AP is Lower


   00 01 02 03 04              00


      11 12 13 14              10 11


         22 23 24              20 21 22


            33 34              30 31 32 33


               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last


  three columns of AP upper. The lower triangle A(3:4,0:1) consists of


  conjugate-transpose of the first two   columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:1,1:2) consists of


  conjugate-transpose of the last two   columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N odd  and TRANSR = 'N'.


         RFP A                   RFP A


                                   -- --


        02 03 04                00 33 43


                                      --


        12 13 14                10 11 44


        22 23 24                20 21 22


        --


        00 33 34                30 31 32


        -- --


        01 11 44                40 41 42


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- --                   -- -- -- -- -- --


     02 12 22 00 01             00 10 20 30 40 50


     -- -- -- --                   -- -- -- -- --


     03 13 23 33 11             33 11 21 31 41 51


     -- -- -- -- --                   -- -- -- --


     04 14 24 34 44             43 44 22 32 42 52

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPFTRI computes the inverse of a complex Hermitian positive definite


 matrix A using the Cholesky factorization A = U**H*U or A = L*L**H


 computed by ZPFTRF.

TRANSR

          TRANSR is CHARACTER*1


          = 'N':  The Normal TRANSR of RFP A is stored;


          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

A

          A is COMPLEX*16 array, dimension ( N*(N+1)/2 );


          On entry, the Hermitian matrix A in RFP format. RFP format is


          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'


          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is


          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'C' then RFP is


          the Conjugate-transpose of RFP A as defined when


          TRANSR = 'N'. The contents of RFP A are defined by UPLO as


          follows: If UPLO = 'U' the RFP A contains the nt elements of


          upper packed A. If UPLO = 'L' the RFP A contains the elements


          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =


          'C'. When TRANSR is 'N' the LDA is N+1 when N is even and N


          is odd. See the Note below for more details.


          On exit, the Hermitian inverse of the original matrix, in the


          same storage format.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the (i,i) element of the factor U or L is


                zero, and the inverse could not be computed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  We first consider Standard Packed Format when N is even.


  We give an example where N = 6.


      AP is Upper             AP is Lower


   00 01 02 03 04 05       00


      11 12 13 14 15       10 11


         22 23 24 25       20 21 22


            33 34 35       30 31 32 33


               44 45       40 41 42 43 44


                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last


  three columns of AP upper. The lower triangle A(4:6,0:2) consists of


  conjugate-transpose of the first three columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:2,0:2) consists of


  conjugate-transpose of the last three columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N even and TRANSR = 'N'.


         RFP A                   RFP A


                                -- -- --


        03 04 05                33 43 53


                                   -- --


        13 14 15                00 44 54


                                      --


        23 24 25                10 11 55


        33 34 35                20 21 22


        --


        00 44 45                30 31 32


        -- --


        01 11 55                40 41 42


        -- -- --


        02 12 22                50 51 52


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- -- --                -- -- -- -- -- --


     03 13 23 33 00 01 02    33 00 10 20 30 40 50


     -- -- -- -- --                -- -- -- -- --


     04 14 24 34 44 11 12    43 44 11 21 31 41 51


     -- -- -- -- -- --                -- -- -- --


     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We next  consider Standard Packed Format when N is odd.


  We give an example where N = 5.


     AP is Upper                 AP is Lower


   00 01 02 03 04              00


      11 12 13 14              10 11


         22 23 24              20 21 22


            33 34              30 31 32 33


               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last


  three columns of AP upper. The lower triangle A(3:4,0:1) consists of


  conjugate-transpose of the first two   columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:1,1:2) consists of


  conjugate-transpose of the last two   columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N odd  and TRANSR = 'N'.


         RFP A                   RFP A


                                   -- --


        02 03 04                00 33 43


                                      --


        12 13 14                10 11 44


        22 23 24                20 21 22


        --


        00 33 34                30 31 32


        -- --


        01 11 44                40 41 42


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- --                   -- -- -- -- -- --


     02 12 22 00 01             00 10 20 30 40 50


     -- -- -- --                   -- -- -- -- --


     03 13 23 33 11             33 11 21 31 41 51


     -- -- -- -- --                   -- -- -- --


     04 14 24 34 44             43 44 22 32 42 52

 ZPFTRS solves a system of linear equations A*X = B with a Hermitian


 positive definite matrix A using the Cholesky factorization


 A = U**H*U or A = L*L**H computed by ZPFTRF.

TRANSR

          TRANSR is CHARACTER*1


          = 'N':  The Normal TRANSR of RFP A is stored;


          = 'C':  The Conjugate-transpose TRANSR of RFP A is stored.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of RFP A is stored;


          = 'L':  Lower triangle of RFP A is stored.

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 matrix B.  NRHS >= 0.

A

          A is COMPLEX*16 array, dimension ( N*(N+1)/2 );


          The triangular factor U or L from the Cholesky factorization


          of RFP A = U**H*U or RFP A = L*L**H, as computed by ZPFTRF.


          See note below for more details about RFP A.

B

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


          On entry, the right hand side matrix B.


          On exit, the solution matrix X.

LDB

          LDB is INTEGER


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

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  We first consider Standard Packed Format when N is even.


  We give an example where N = 6.


      AP is Upper             AP is Lower


   00 01 02 03 04 05       00


      11 12 13 14 15       10 11


         22 23 24 25       20 21 22


            33 34 35       30 31 32 33


               44 45       40 41 42 43 44


                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last


  three columns of AP upper. The lower triangle A(4:6,0:2) consists of


  conjugate-transpose of the first three columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:2,0:2) consists of


  conjugate-transpose of the last three columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N even and TRANSR = 'N'.


         RFP A                   RFP A


                                -- -- --


        03 04 05                33 43 53


                                   -- --


        13 14 15                00 44 54


                                      --


        23 24 25                10 11 55


        33 34 35                20 21 22


        --


        00 44 45                30 31 32


        -- --


        01 11 55                40 41 42


        -- -- --


        02 12 22                50 51 52


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- -- --                -- -- -- -- -- --


     03 13 23 33 00 01 02    33 00 10 20 30 40 50


     -- -- -- -- --                -- -- -- -- --


     04 14 24 34 44 11 12    43 44 11 21 31 41 51


     -- -- -- -- -- --                -- -- -- --


     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We next  consider Standard Packed Format when N is odd.


  We give an example where N = 5.


     AP is Upper                 AP is Lower


   00 01 02 03 04              00


      11 12 13 14              10 11


         22 23 24              20 21 22


            33 34              30 31 32 33


               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:


  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last


  three columns of AP upper. The lower triangle A(3:4,0:1) consists of


  conjugate-transpose of the first two   columns of AP upper.


  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first


  three columns of AP lower. The upper triangle A(0:1,1:2) consists of


  conjugate-transpose of the last two   columns of AP lower.


  To denote conjugate we place -- above the element. This covers the


  case N odd  and TRANSR = 'N'.


         RFP A                   RFP A


                                   -- --


        02 03 04                00 33 43


                                      --


        12 13 14                10 11 44


        22 23 24                20 21 22


        --


        00 33 34                30 31 32


        -- --


        01 11 44                40 41 42


  Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate-


  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A


     -- -- --                   -- -- -- -- -- --


     02 12 22 00 01             00 10 20 30 40 50


     -- -- -- --                   -- -- -- -- --


     03 13 23 33 11             33 11 21 31 41 51


     -- -- -- -- --                   -- -- -- --


     04 14 24 34 44             43 44 22 32 42 52

 ZPPCON estimates the reciprocal of the condition number (in the


 1-norm) of a complex Hermitian positive definite packed matrix using


 the Cholesky factorization A = U**H*U or A = L*L**H computed by


 ZPPTRF.


 An estimate is obtained for norm(inv(A)), and the reciprocal of the


 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The triangular factor U or L from the Cholesky factorization


          A = U**H*U or A = L*L**H, packed columnwise in a linear


          array.  The j-th column of U or L is stored in the array AP


          as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.

ANORM

          ANORM is DOUBLE PRECISION


          The 1-norm (or infinity-norm) of the Hermitian matrix A.

RCOND

          RCOND is DOUBLE PRECISION


          The reciprocal of the condition number of the matrix A,


          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an


          estimate of the 1-norm of inv(A) computed in this routine.

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPPEQU computes row and column scalings intended to equilibrate a


 Hermitian positive definite matrix A in packed storage and reduce


 its condition number (with respect to the two-norm).  S contains the


 scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix


 B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal.


 This choice of S puts the condition number of B within a factor N of


 the smallest possible condition number over all possible diagonal


 scalings.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The upper or lower 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.

S

          S is DOUBLE PRECISION array, dimension (N)


          If INFO = 0, S contains the scale factors for A.

SCOND

          SCOND is DOUBLE PRECISION


          If INFO = 0, S contains the ratio of the smallest S(i) to


          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too


          large nor too small, it is not worth scaling by S.

AMAX

          AMAX is DOUBLE PRECISION


          Absolute value of largest matrix element.  If AMAX is very


          close to overflow or very close to underflow, the matrix


          should be scaled.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the i-th diagonal element is nonpositive.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPPRFS improves the computed solution to a system of linear


 equations when the coefficient matrix is Hermitian positive definite


 and packed, and provides error bounds and backward error estimates


 for the solution.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

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 B and X.  NRHS >= 0.

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The upper or lower triangle of the Hermitian matrix A, packed


          columnwise in a linear array.  The j-th column of A is stored


          in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

AFP

          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The triangular factor U or L from the Cholesky factorization


          A = U**H*U or A = L*L**H, as computed by DPPTRF/ZPPTRF,


          packed columnwise in a linear array in the same format as A


          (see AP).

B

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


          The right hand side matrix B.

LDB

          LDB is INTEGER


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

X

          X is COMPLEX*16 array, dimension (LDX,NRHS)


          On entry, the solution matrix X, as computed by ZPPTRS.


          On exit, the improved solution matrix X.

LDX

          LDX is INTEGER


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

FERR

          FERR is DOUBLE PRECISION array, dimension (NRHS)


          The estimated forward error bound for each solution vector


          X(j) (the j-th column of the solution matrix X).


          If XTRUE is the true solution corresponding to X(j), FERR(j)


          is an estimated upper bound for the magnitude of the largest


          element in (X(j) - XTRUE) divided by the magnitude of the


          largest element in X(j).  The estimate is as reliable as


          the estimate for RCOND, and is almost always a slight


          overestimate of the true error.

BERR

          BERR is DOUBLE PRECISION array, dimension (NRHS)


          The componentwise relative backward error of each solution


          vector X(j) (i.e., the smallest relative change in


          any element of A or B that makes X(j) an exact solution).

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

  ITMAX is the maximum number of steps of iterative refinement.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPPTRF computes the Cholesky factorization of a complex Hermitian


 positive definite matrix A stored in packed format.


 The factorization has the form


    A = U**H * U,  if UPLO = 'U', or


    A = L  * L**H,  if UPLO = 'L',


 where U is an upper triangular matrix and L is lower triangular.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the upper or lower triangle of the Hermitian matrix


          A, packed columnwise in a linear array.  The j-th column of A


          is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.


          See below for further details.


          On exit, if INFO = 0, the triangular factor U or L from the


          Cholesky factorization A = U**H*U or A = L*L**H, in the same


          storage format as A.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the leading minor of order i is not


                positive definite, and the factorization could not be


                completed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The packed storage scheme is illustrated by the following example


  when N = 4, UPLO = 'U':


  Two-dimensional storage of the Hermitian matrix A:


     a11 a12 a13 a14


         a22 a23 a24


             a33 a34     (aij = conjg(aji))


                 a44


  Packed storage of the upper triangle of A:


  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]

 ZPPTRI computes the inverse of a complex Hermitian positive definite


 matrix A using the Cholesky factorization A = U**H*U or A = L*L**H


 computed by ZPPTRF.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangular factor is stored in AP;


          = 'L':  Lower triangular factor is stored in AP.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the triangular factor U or L from the Cholesky


          factorization A = U**H*U or A = L*L**H, packed columnwise as


          a linear array.  The j-th column of U or L is stored in the


          array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.


          On exit, the upper or lower triangle of the (Hermitian)


          inverse of A, overwriting the input factor U or L.

INFO

          INFO is INTEGER


          = 0:  successful exit


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


          > 0:  if INFO = i, the (i,i) element of the factor U or L is


                zero, and the inverse could not be computed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPPTRS solves a system of linear equations A*X = B with a Hermitian


 positive definite matrix A in packed storage using the Cholesky


 factorization A = U**H * U or A = L * L**H computed by ZPPTRF.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

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 matrix B.  NRHS >= 0.

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The triangular factor U or L from the Cholesky factorization


          A = U**H * U or A = L * L**H, packed columnwise in a linear


          array.  The j-th column of U or L is stored in the array AP


          as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.

B

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


          On entry, the right hand side matrix B.


          On exit, the solution matrix X.

LDB

          LDB is INTEGER


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

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPSTF2 computes the Cholesky factorization with complete


 pivoting of a complex Hermitian positive semidefinite matrix A.


 The factorization has the form


    P**T * A * P = U**H * U ,  if UPLO = 'U',


    P**T * A * P = L  * L**H,  if UPLO = 'L',


 where U is an upper triangular matrix and L is lower triangular, and


 P is stored as vector PIV.


 This algorithm does not attempt to check that A is positive


 semidefinite. This version of the algorithm calls level 2 BLAS.

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.

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, the symmetric matrix A.  If UPLO = 'U', the leading


          n by n upper triangular part of A contains the upper


          triangular part of the matrix A, and the strictly lower


          triangular part of A is not referenced.  If UPLO = 'L', the


          leading n by n lower triangular part of A contains the lower


          triangular part of the matrix A, and the strictly upper


          triangular part of A is not referenced.


          On exit, if INFO = 0, the factor U or L from the Cholesky


          factorization as above.

PIV

          PIV is INTEGER array, dimension (N)


          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.

RANK

          RANK is INTEGER


          The rank of A given by the number of steps the algorithm


          completed.

TOL

          TOL is DOUBLE PRECISION


          User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )


          will be used. The algorithm terminates at the (K-1)st step


          if the pivot <= TOL.

LDA

          LDA is INTEGER


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

WORK

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


          Work space.

INFO

          INFO is INTEGER


          < 0: If INFO = -K, the K-th argument had an illegal value,


          = 0: algorithm completed successfully, and


          > 0: the matrix A is either rank deficient with computed rank


               as returned in RANK, or is not positive semidefinite. See


               Section 7 of LAPACK Working Note #161 for further


               information.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZPSTRF computes the Cholesky factorization with complete


 pivoting of a complex Hermitian positive semidefinite matrix A.


 The factorization has the form


    P**T * A * P = U**H * U ,  if UPLO = 'U',


    P**T * A * P = L  * L**H,  if UPLO = 'L',


 where U is an upper triangular matrix and L is lower triangular, and


 P is stored as vector PIV.


 This algorithm does not attempt to check that A is positive


 semidefinite. This version of the algorithm calls level 3 BLAS.

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.

A

          A is COMPLEX*16 array, dimension (LDA,N)


          On entry, the symmetric matrix A.  If UPLO = 'U', the leading


          n by n upper triangular part of A contains the upper


          triangular part of the matrix A, and the strictly lower


          triangular part of A is not referenced.  If UPLO = 'L', the


          leading n by n lower triangular part of A contains the lower


          triangular part of the matrix A, and the strictly upper


          triangular part of A is not referenced.


          On exit, if INFO = 0, the factor U or L from the Cholesky


          factorization as above.

LDA

          LDA is INTEGER


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

PIV

          PIV is INTEGER array, dimension (N)


          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.

RANK

          RANK is INTEGER


          The rank of A given by the number of steps the algorithm


          completed.

TOL

          TOL is DOUBLE PRECISION


          User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )


          will be used. The algorithm terminates at the (K-1)st step


          if the pivot <= TOL.

WORK

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


          Work space.

INFO

          INFO is INTEGER


          < 0: If INFO = -K, the K-th argument had an illegal value,


          = 0: algorithm completed successfully, and


          > 0: the matrix A is either rank deficient with computed rank


               as returned in RANK, or is not positive semidefinite. See


               Section 7 of LAPACK Working Note #161 for further


               information.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZSPCON estimates the reciprocal of the condition number (in the


 1-norm) of a complex symmetric packed matrix A using the


 factorization A = U*D*U**T or A = L*D*L**T computed by ZSPTRF.


 An estimate is obtained for norm(inv(A)), and the reciprocal of the


 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

UPLO

          UPLO is CHARACTER*1


          Specifies whether the details of the factorization are stored


          as an upper or lower triangular matrix.


          = 'U':  Upper triangular, form is A = U*D*U**T;


          = 'L':  Lower triangular, form is A = L*D*L**T.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The block diagonal matrix D and the multipliers used to


          obtain the factor U or L as computed by ZSPTRF, stored as a


          packed triangular matrix.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZSPTRF.

ANORM

          ANORM is DOUBLE PRECISION


          The 1-norm of the original matrix A.

RCOND

          RCOND is DOUBLE PRECISION


          The reciprocal of the condition number of the matrix A,


          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an


          estimate of the 1-norm of inv(A) computed in this routine.

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZSPRFS improves the computed solution to a system of linear


 equations when the coefficient matrix is symmetric indefinite


 and packed, and provides error bounds and backward error estimates


 for the solution.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

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 B and X.  NRHS >= 0.

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The upper or lower triangle of the symmetric matrix A, packed


          columnwise in a linear array.  The j-th column of A is stored


          in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

AFP

          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)


          The factored form of the matrix A.  AFP contains the block


          diagonal matrix D and the multipliers used to obtain the


          factor U or L from the factorization A = U*D*U**T or


          A = L*D*L**T as computed by ZSPTRF, stored as a packed


          triangular matrix.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZSPTRF.

B

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


          The right hand side matrix B.

LDB

          LDB is INTEGER


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

X

          X is COMPLEX*16 array, dimension (LDX,NRHS)


          On entry, the solution matrix X, as computed by ZSPTRS.


          On exit, the improved solution matrix X.

LDX

          LDX is INTEGER


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

FERR

          FERR is DOUBLE PRECISION array, dimension (NRHS)


          The estimated forward error bound for each solution vector


          X(j) (the j-th column of the solution matrix X).


          If XTRUE is the true solution corresponding to X(j), FERR(j)


          is an estimated upper bound for the magnitude of the largest


          element in (X(j) - XTRUE) divided by the magnitude of the


          largest element in X(j).  The estimate is as reliable as


          the estimate for RCOND, and is almost always a slight


          overestimate of the true error.

BERR

          BERR is DOUBLE PRECISION array, dimension (NRHS)


          The componentwise relative backward error of each solution


          vector X(j) (i.e., the smallest relative change in


          any element of A or B that makes X(j) an exact solution).

WORK

          WORK is COMPLEX*16 array, dimension (2*N)

RWORK

          RWORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

  ITMAX is the maximum number of steps of iterative refinement.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 ZSPTRF computes the factorization of a complex symmetric matrix A


 stored in packed format using the Bunch-Kaufman diagonal pivoting


 method:


    A = U*D*U**T  or  A = L*D*L**T


 where U (or L) is a product of permutation and unit upper (lower)


 triangular matrices, and D is symmetric and block diagonal with


 1-by-1 and 2-by-2 diagonal blocks.

UPLO

          UPLO is CHARACTER*1


          = 'U':  Upper triangle of A is stored;


          = 'L':  Lower triangle of A is stored.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the upper or lower triangle of the symmetric matrix


          A, packed columnwise in a linear array.  The j-th column of A


          is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;


          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.


          On exit, the block diagonal matrix D and the multipliers used


          to obtain the factor U or L, stored as a packed triangular


          matrix overwriting A (see below for further details).

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D.


          If IPIV(k) > 0, then rows and columns k and IPIV(k) were


          interchanged and D(k,k) is a 1-by-1 diagonal block.


          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and


          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)


          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =


          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were


          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.

INFO

          INFO is INTEGER


          = 0: successful exit


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


          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization


               has been completed, but the block diagonal matrix D is


               exactly singular, and division by zero will occur if it


               is used to solve a system of equations.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  5-96 - Based on modifications by J. Lewis, Boeing Computer Services


         Company


  If UPLO = 'U', then A = U*D*U**T, where


     U = P(n)*U(n)* ... *P(k)U(k)* ...,


  i.e., U is a product of terms P(k)*U(k), where k decreases from n to


  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1


  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as


  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such


  that if the diagonal block D(k) is of order s (s = 1 or 2), then


             (   I    v    0   )   k-s


     U(k) =  (   0    I    0   )   s


             (   0    0    I   )   n-k


                k-s   s   n-k


  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).


  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),


  and A(k,k), and v overwrites A(1:k-2,k-1:k).


  If UPLO = 'L', then A = L*D*L**T, where


     L = P(1)*L(1)* ... *P(k)*L(k)* ...,


  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to


  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1


  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as


  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such


  that if the diagonal block D(k) is of order s (s = 1 or 2), then


             (   I    0     0   )  k-1


     L(k) =  (   0    I     0   )  s


             (   0    v     I   )  n-k-s+1


                k-1   s  n-k-s+1


  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).


  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),


  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).

 ZSPTRI computes the inverse of a complex symmetric indefinite matrix


 A in packed storage using the factorization A = U*D*U**T or


 A = L*D*L**T computed by ZSPTRF.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the details of the factorization are stored


          as an upper or lower triangular matrix.


          = 'U':  Upper triangular, form is A = U*D*U**T;


          = 'L':  Lower triangular, form is A = L*D*L**T.

N

          N is INTEGER


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

AP

          AP is COMPLEX*16 array, dimension (N*(N+1)/2)


          On entry, the block diagonal matrix D and the multipliers


          used to obtain the factor U or L as computed by ZSPTRF,


          stored as a packed triangular matrix.


          On exit, if INFO = 0, the (symmetric) inverse of the original


          matrix, stored as a packed triangular matrix. The j-th column


          of inv(A) is stored in the array AP as follows:


          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;


          if UPLO = 'L',


             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.

IPIV

          IPIV is INTEGER array, dimension (N)


          Details of the interchanges and the block structure of D


          as determined by ZSPTRF.

WORK

          WORK is COMPLEX*16 array, dimension (N)

INFO

          INFO is INTEGER


          = 0: successful exit


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


          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its


               inverse could not be computed.