Provided by: liblapack-doc_3.8.0-2_all

**NAME**

complex16OTHERauxiliary

**SYNOPSIS**

Functionssubroutineclag2z(M, N, SA, LDSA, A, LDA, INFO)CLAG2Zconverts a complex single precision matrix to a complex double precision matrix. double precision functiondzsum1(N, CX, INCX)DZSUM1forms the 1-norm of the complex vector using the true absolute value. integer functionilazlc(M, N, A, LDA)ILAZLCscans a matrix for its last non-zero column. integer functionilazlr(M, N, A, LDA)ILAZLRscans a matrix for its last non-zero row. subroutinezdrscl(N, SA, SX, INCX)ZDRSCLmultiplies a vector by the reciprocal of a real scalar. subroutinezlabrd(M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)ZLABRDreduces the first nb rows and columns of a general matrix to a bidiagonal form. subroutinezlacgv(N, X, INCX)ZLACGVconjugates a complex vector. subroutinezlacn2(N, V, X, EST, KASE, ISAVE)ZLACN2estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. subroutinezlacon(N, V, X, EST, KASE)ZLACONestimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products. subroutinezlacp2(UPLO, M, N, A, LDA, B, LDB)ZLACP2copies all or part of a real two-dimensional array to a complex array. subroutinezlacpy(UPLO, M, N, A, LDA, B, LDB)ZLACPYcopies all or part of one two-dimensional array to another. subroutinezlacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)ZLACRMmultiplies a complex matrix by a square real matrix. subroutinezlacrt(N, CX, INCX, CY, INCY, C, S)ZLACRTperforms a linear transformation of a pair of complex vectors. complex *16 functionzladiv(X, Y)ZLADIVperforms complex division in real arithmetic, avoiding unnecessary overflow. subroutinezlaein(RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, EPS3, SMLNUM, INFO)ZLAEINcomputes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. subroutinezlaev2(A, B, C, RT1, RT2, CS1, SN1)ZLAEV2computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix. subroutinezlag2c(M, N, A, LDA, SA, LDSA, INFO)ZLAG2Cconverts a complex double precision matrix to a complex single precision matrix. subroutinezlags2(UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)ZLAGS2subroutinezlagtm(TRANS, N, NRHS, ALPHA, DL, D, DU, X, LDX, BETA, B, LDB)ZLAGTMperforms a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix, B and C are rectangular matrices, and α and β are scalars, which may be 0, 1, or -1. subroutinezlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)ZLAHQRcomputes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm. subroutinezlahr2(N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)ZLAHR2reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A. subroutinezlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)ZLAIC1applies one step of incremental condition estimation. double precision functionzlangt(NORM, N, DL, D, DU)ZLANGTreturns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix. double precision functionzlanhb(NORM, UPLO, N, K, AB, LDAB, WORK)ZLANHBreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix. double precision functionzlanhp(NORM, UPLO, N, AP, WORK)ZLANHPreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form. double precision functionzlanhs(NORM, N, A, LDA, WORK)ZLANHSreturns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix. double precision functionzlanht(NORM, N, D, E)ZLANHTreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix. double precision functionzlansb(NORM, UPLO, N, K, AB, LDAB, WORK)ZLANSBreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix. double precision functionzlansp(NORM, UPLO, N, AP, WORK)ZLANSPreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form. double precision functionzlantb(NORM, UPLO, DIAG, N, K, AB, LDAB, WORK)ZLANTBreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular band matrix. double precision functionzlantp(NORM, UPLO, DIAG, N, AP, WORK)ZLANTPreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form. double precision functionzlantr(NORM, UPLO, DIAG, M, N, A, LDA, WORK)ZLANTRreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix. subroutinezlapll(N, X, INCX, Y, INCY, SSMIN)ZLAPLLmeasures the linear dependence of two vectors. subroutinezlapmr(FORWRD, M, N, X, LDX, K)ZLAPMRrearranges rows of a matrix as specified by a permutation vector. subroutinezlapmt(FORWRD, M, N, X, LDX, K)ZLAPMTperforms a forward or backward permutation of the columns of a matrix. subroutinezlaqhb(UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)ZLAQHBscales a Hermitian band matrix, using scaling factors computed by cpbequ. subroutinezlaqhp(UPLO, N, AP, S, SCOND, AMAX, EQUED)ZLAQHPscales a Hermitian matrix stored in packed form. subroutinezlaqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)ZLAQP2computes a QR factorization with column pivoting of the matrix block. subroutinezlaqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)ZLAQPScomputes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3. subroutinezlaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)ZLAQR0computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. subroutinezlaqr1(N, H, LDH, S1, S2, V)ZLAQR1sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts. subroutinezlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)ZLAQR2performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). subroutinezlaqr3(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)ZLAQR3performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). subroutinezlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)ZLAQR4computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition. subroutinezlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)ZLAQR5performs a single small-bulge multi-shift QR sweep. subroutinezlaqsb(UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED)ZLAQSBscales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ. subroutinezlaqsp(UPLO, N, AP, S, SCOND, AMAX, EQUED)ZLAQSPscales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ. subroutinezlar1v(N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)ZLAR1Vcomputes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI. subroutinezlar2v(N, X, Y, Z, INCX, C, S, INCC)ZLAR2Vapplies a vector of plane rotations with real cosines and complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices. subroutinezlarcm(M, N, A, LDA, B, LDB, C, LDC, RWORK)ZLARCMcopies all or part of a real two-dimensional array to a complex array. subroutinezlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)ZLARFapplies an elementary reflector to a general rectangular matrix. subroutinezlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)ZLARFBapplies a block reflector or its conjugate-transpose to a general rectangular matrix. subroutinezlarfg(N, ALPHA, X, INCX, TAU)ZLARFGgenerates an elementary reflector (Householder matrix). subroutinezlarfgp(N, ALPHA, X, INCX, TAU)ZLARFGPgenerates an elementary reflector (Householder matrix) with non-negative beta. subroutinezlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)ZLARFTforms the triangular factor T of a block reflector H = I - vtvH subroutinezlarfx(SIDE, M, N, V, TAU, C, LDC, WORK)ZLARFXapplies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10. subroutinezlargv(N, X, INCX, Y, INCY, C, INCC)ZLARGVgenerates a vector of plane rotations with real cosines and complex sines. subroutinezlarnv(IDIST, ISEED, N, X)ZLARNVreturns a vector of random numbers from a uniform or normal distribution. subroutinezlarrv(N, VL, VU, D, L, PIVMIN, ISPLIT, M, DOL, DOU, MINRGP, RTOL1, RTOL2, W, WERR, WGAP, IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, WORK, IWORK, INFO)ZLARRVcomputes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. subroutinezlartg(F, G, CS, SN, R)ZLARTGgenerates a plane rotation with real cosine and complex sine. subroutinezlartv(N, X, INCX, Y, INCY, C, S, INCC)ZLARTVapplies a vector of plane rotations with real cosines and complex sines to the elements of a pair of vectors. subroutinezlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)ZLASCLmultiplies a general rectangular matrix by a real scalar defined as cto/cfrom. subroutinezlaset(UPLO, M, N, ALPHA, BETA, A, LDA)ZLASETinitializes the off-diagonal elements and the diagonal elements of a matrix to given values. subroutinezlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)ZLASRapplies a sequence of plane rotations to a general rectangular matrix. subroutinezlassq(N, X, INCX, SCALE, SUMSQ)ZLASSQupdates a sum of squares represented in scaled form. subroutinezlaswp(N, A, LDA, K1, K2, IPIV, INCX)ZLASWPperforms a series of row interchanges on a general rectangular matrix. subroutinezlat2c(UPLO, N, A, LDA, SA, LDSA, INFO)ZLAT2Cconverts a double complex triangular matrix to a complex triangular matrix. subroutinezlatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)ZLATBSsolves a triangular banded system of equations. subroutinezlatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)ZLATDFuses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate. subroutinezlatps(UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE, CNORM, INFO)ZLATPSsolves a triangular system of equations with the matrix held in packed storage. subroutinezlatrd(UPLO, N, NB, A, LDA, E, TAU, W, LDW)ZLATRDreduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an unitary similarity transformation. subroutinezlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)ZLATRSsolves a triangular system of equations with the scale factor set to prevent overflow. subroutinezlauu2(UPLO, N, A, LDA, INFO)ZLAUU2computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm). subroutinezlauum(UPLO, N, A, LDA, INFO)ZLAUUMcomputes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). subroutinezrot(N, CX, INCX, CY, INCY, C, S)ZROTapplies a plane rotation with real cosine and complex sine to a pair of complex vectors. subroutinezspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)ZSPMVcomputes a matrix-vector product for complex vectors using a complex symmetric packed matrix subroutinezspr(UPLO, N, ALPHA, X, INCX, AP)ZSPRperforms the symmetrical rank-1 update of a complex symmetric packed matrix. subroutineztprfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK)ZTPRFBapplies a real or complex 'triangular-pentagonal' blocked reflector to a real or complex matrix, which is composed of two blocks.

**Detailed** **Description**

This is the group of complex16 other auxiliary routines

**Function** **Documentation**

subroutineclag2z(integerM,integerN,complex,dimension(ldsa,*)SA,integerLDSA,complex*16,dimension(lda,*)A,integerLDA,integerINFO)CLAG2Zconverts a complex single precision matrix to a complex double precision matrix.Purpose:CLAG2Z converts a COMPLEX matrix, SA, to a COMPLEX*16 matrix, A. Note that while it is possible to overflow while converting from double to single, it is not possible to overflow when converting from single to double. This is an auxiliary routine so there is no argument checking.Parameters:MM is INTEGER The number of lines of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0.SASA is COMPLEX array, dimension (LDSA,N) On entry, the M-by-N coefficient matrix SA.LDSALDSA is INTEGER The leading dimension of the array SA. LDSA >= max(1,M).AA is COMPLEX*16 array, dimension (LDA,N) On exit, the M-by-N coefficient matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).INFOINFO is INTEGER = 0: successful exitAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctiondzsum1(integerN,complex*16,dimension(*)CX,integerINCX)DZSUM1forms the 1-norm of the complex vector using the true absolute value.Purpose:DZSUM1 takes the sum of the absolute values of a complex vector and returns a double precision result. Based on DZASUM from the Level 1 BLAS. The change is to use the 'genuine' absolute value.Parameters:NN is INTEGER The number of elements in the vector CX.CXCX is COMPLEX*16 array, dimension (N) The vector whose elements will be summed.INCXINCX is INTEGER The spacing between successive values of CX. INCX > 0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Nick Higham for use with ZLACON.integerfunctionilazlc(integerM,integerN,complex*16,dimension(lda,*)A,integerLDA)ILAZLCscans a matrix for its last non-zero column.Purpose:ILAZLC scans A for its last non-zero column.Parameters:MM is INTEGER The number of rows of the matrix A.NN is INTEGER The number of columns of the matrix A.AA is COMPLEX*16 array, dimension (LDA,N) The m by n matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016integerfunctionilazlr(integerM,integerN,complex*16,dimension(lda,*)A,integerLDA)ILAZLRscans a matrix for its last non-zero row.Purpose:ILAZLR scans A for its last non-zero row.Parameters:MM is INTEGER The number of rows of the matrix A.NN is INTEGER The number of columns of the matrix A.AA is COMPLEX*16 array, dimension (LDA,N) The m by n matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezdrscl(integerN,doubleprecisionSA,complex*16,dimension(*)SX,integerINCX)ZDRSCLmultiplies a vector by the reciprocal of a real scalar.Purpose:ZDRSCL multiplies an n-element complex vector x by the real scalar 1/a. This is done without overflow or underflow as long as the final result x/a does not overflow or underflow.Parameters:NN is INTEGER The number of components of the vector x.SASA is DOUBLE PRECISION The scalar a which is used to divide each component of x. SA must be >= 0, or the subroutine will divide by zero.SXSX is COMPLEX*16 array, dimension (1+(N-1)*abs(INCX)) The n-element vector x.INCXINCX is INTEGER The increment between successive values of the vector SX. > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= nAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlabrd(integerM,integerN,integerNB,complex*16,dimension(lda,*)A,integerLDA,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,complex*16,dimension(*)TAUQ,complex*16,dimension(*)TAUP,complex*16,dimension(ldx,*)X,integerLDX,complex*16,dimension(ldy,*)Y,integerLDY)ZLABRDreduces the first nb rows and columns of a general matrix to a bidiagonal form.Purpose:ZLABRD reduces the first NB rows and columns of a complex general m by n matrix A to upper or lower real bidiagonal form by a unitary transformation Q**H * A * P, and returns the matrices X and Y which are needed to apply the transformation to the unreduced part of A. If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower bidiagonal form. This is an auxiliary routine called by ZGEBRDParameters:MM is INTEGER The number of rows in the matrix A.NN is INTEGER The number of columns in the matrix A.NBNB is INTEGER The number of leading rows and columns of A to be reduced.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the m by n general matrix to be reduced. On exit, the first NB rows and columns of the matrix are overwritten; the rest of the array is unchanged. If m >= n, elements on and below the diagonal in the first NB columns, with the array TAUQ, represent the unitary matrix Q as a product of elementary reflectors; and elements above the diagonal in the first NB rows, with the array TAUP, represent the unitary matrix P as a product of elementary reflectors. If m < n, elements below the diagonal in the first NB columns, with the array TAUQ, represent the unitary matrix Q as a product of elementary reflectors, and elements on and above the diagonal in the first NB rows, with the array TAUP, represent the unitary matrix P as a product of elementary reflectors. See Further Details.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).DD is DOUBLE PRECISION array, dimension (NB) The diagonal elements of the first NB rows and columns of the reduced matrix. D(i) = A(i,i).EE is DOUBLE PRECISION array, dimension (NB) The off-diagonal elements of the first NB rows and columns of the reduced matrix.TAUQTAUQ is COMPLEX*16 array, dimension (NB) The scalar factors of the elementary reflectors which represent the unitary matrix Q. See Further Details.TAUPTAUP is COMPLEX*16 array, dimension (NB) The scalar factors of the elementary reflectors which represent the unitary matrix P. See Further Details.XX is COMPLEX*16 array, dimension (LDX,NB) The m-by-nb matrix X required to update the unreduced part of A.LDXLDX is INTEGER The leading dimension of the array X. LDX >= max(1,M).YY is COMPLEX*16 array, dimension (LDY,NB) The n-by-nb matrix Y required to update the unreduced part of A.LDYLDY is INTEGER The leading dimension of the array Y. LDY >= max(1,N).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017FurtherDetails:The matrices Q and P are represented as products of elementary reflectors: Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) Each H(i) and G(i) has the form: H(i) = I - tauq * v * v**H and G(i) = I - taup * u * u**H where tauq and taup are complex scalars, and v and u are complex vectors. If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). The elements of the vectors v and u together form the m-by-nb matrix V and the nb-by-n matrix U**H which are needed, with X and Y, to apply the transformation to the unreduced part of the matrix, using a block update of the form: A := A - V*Y**H - X*U**H. The contents of A on exit are illustrated by the following examples with nb = 2: m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) ( v1 v2 a a a ) ( v1 1 a a a a ) ( v1 v2 a a a ) ( v1 v2 a a a a ) ( v1 v2 a a a ) ( v1 v2 a a a a ) ( v1 v2 a a a ) where a denotes an element of the original matrix which is unchanged, vi denotes an element of the vector defining H(i), and ui an element of the vector defining G(i).subroutinezlacgv(integerN,complex*16,dimension(*)X,integerINCX)ZLACGVconjugates a complex vector.Purpose:ZLACGV conjugates a complex vector of length N.Parameters:NN is INTEGER The length of the vector X. N >= 0.XX is COMPLEX*16 array, dimension (1+(N-1)*abs(INCX)) On entry, the vector of length N to be conjugated. On exit, X is overwritten with conjg(X).INCXINCX is INTEGER The spacing between successive elements of X.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlacn2(integerN,complex*16,dimension(*)V,complex*16,dimension(*)X,doubleprecisionEST,integerKASE,integer,dimension(3)ISAVE)ZLACN2estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products.Purpose:ZLACN2 estimates the 1-norm of a square, complex matrix A. Reverse communication is used for evaluating matrix-vector products.Parameters:NN is INTEGER The order of the matrix. N >= 1.VV is COMPLEX*16 array, dimension (N) On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned).XX is COMPLEX*16 array, dimension (N) On an intermediate return, X should be overwritten by A * X, if KASE=1, A**H * X, if KASE=2, where A**H is the conjugate transpose of A, and ZLACN2 must be re-called with all the other parameters unchanged.ESTEST is DOUBLE PRECISION On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be unchanged from the previous call to ZLACN2. On exit, EST is an estimate (a lower bound) for norm(A).KASEKASE is INTEGER On the initial call to ZLACN2, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether X should be overwritten by A * X or A**H * X. On the final return from ZLACN2, KASE will again be 0.ISAVEISAVE is INTEGER array, dimension (3) ISAVE is used to save variables between calls to ZLACN2Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:Originally named CONEST, dated March 16, 1988. Last modified: April, 1999 This is a thread safe version of ZLACON, which uses the array ISAVE in place of a SAVE statement, as follows: ZLACON ZLACN2 JUMP ISAVE(1) J ISAVE(2) ITER ISAVE(3)Contributors:Nick Higham, University of ManchesterReferences:N.J. Higham, 'FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.subroutinezlacon(integerN,complex*16,dimension(n)V,complex*16,dimension(n)X,doubleprecisionEST,integerKASE)ZLACONestimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products.Purpose:ZLACON estimates the 1-norm of a square, complex matrix A. Reverse communication is used for evaluating matrix-vector products.Parameters:NN is INTEGER The order of the matrix. N >= 1.VV is COMPLEX*16 array, dimension (N) On the final return, V = A*W, where EST = norm(V)/norm(W) (W is not returned).XX is COMPLEX*16 array, dimension (N) On an intermediate return, X should be overwritten by A * X, if KASE=1, A**H * X, if KASE=2, where A**H is the conjugate transpose of A, and ZLACON must be re-called with all the other parameters unchanged.ESTEST is DOUBLE PRECISION On entry with KASE = 1 or 2 and JUMP = 3, EST should be unchanged from the previous call to ZLACON. On exit, EST is an estimate (a lower bound) for norm(A).KASEKASE is INTEGER On the initial call to ZLACON, KASE should be 0. On an intermediate return, KASE will be 1 or 2, indicating whether X should be overwritten by A * X or A**H * X. On the final return from ZLACON, KASE will again be 0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:Originally named CONEST, dated March 16, 1988. Last modified: April, 1999Contributors:Nick Higham, University of ManchesterReferences:N.J. Higham, 'FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation', ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.subroutinezlacp2(characterUPLO,integerM,integerN,doubleprecision,dimension(lda,*)A,integerLDA,complex*16,dimension(ldb,*)B,integerLDB)ZLACP2copies all or part of a real two-dimensional array to a complex array.Purpose:ZLACP2 copies all or part of a real two-dimensional matrix A to a complex matrix B.Parameters:UPLOUPLO is CHARACTER*1 Specifies the part of the matrix A to be copied to B. = 'U': Upper triangular part = 'L': Lower triangular part Otherwise: All of the matrix AMM is INTEGER The number of rows of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0.AA is DOUBLE PRECISION array, dimension (LDA,N) The m by n matrix A. If UPLO = 'U', only the upper trapezium is accessed; if UPLO = 'L', only the lower trapezium is accessed.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).BB is COMPLEX*16 array, dimension (LDB,N) On exit, B = A in the locations specified by UPLO.LDBLDB is INTEGER The leading dimension of the array B. LDB >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlacpy(characterUPLO,integerM,integerN,complex*16,dimension(lda,*)A,integerLDA,complex*16,dimension(ldb,*)B,integerLDB)ZLACPYcopies all or part of one two-dimensional array to another.Purpose:ZLACPY copies all or part of a two-dimensional matrix A to another matrix B.Parameters:UPLOUPLO is CHARACTER*1 Specifies the part of the matrix A to be copied to B. = 'U': Upper triangular part = 'L': Lower triangular part Otherwise: All of the matrix AMM is INTEGER The number of rows of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) The m by n matrix A. If UPLO = 'U', only the upper trapezium is accessed; if UPLO = 'L', only the lower trapezium is accessed.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).BB is COMPLEX*16 array, dimension (LDB,N) On exit, B = A in the locations specified by UPLO.LDBLDB is INTEGER The leading dimension of the array B. LDB >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlacrm(integerM,integerN,complex*16,dimension(lda,*)A,integerLDA,doubleprecision,dimension(ldb,*)B,integerLDB,complex*16,dimension(ldc,*)C,integerLDC,doubleprecision,dimension(*)RWORK)ZLACRMmultiplies a complex matrix by a square real matrix.Purpose:ZLACRM performs a very simple matrix-matrix multiplication: C := A * B, where A is M by N and complex; B is N by N and real; C is M by N and complex.Parameters:MM is INTEGER The number of rows of the matrix A and of the matrix C. M >= 0.NN is INTEGER The number of columns and rows of the matrix B and the number of columns of the matrix C. N >= 0.AA is COMPLEX*16 array, dimension (LDA, N) On entry, A contains the M by N matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >=max(1,M).BB is DOUBLE PRECISION array, dimension (LDB, N) On entry, B contains the N by N matrix B.LDBLDB is INTEGER The leading dimension of the array B. LDB >=max(1,N).CC is COMPLEX*16 array, dimension (LDC, N) On exit, C contains the M by N matrix C.LDCLDC is INTEGER The leading dimension of the array C. LDC >=max(1,N).RWORKRWORK is DOUBLE PRECISION array, dimension (2*M*N)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlacrt(integerN,complex*16,dimension(*)CX,integerINCX,complex*16,dimension(*)CY,integerINCY,complex*16C,complex*16S)ZLACRTperforms a linear transformation of a pair of complex vectors.Purpose:ZLACRT performs the operation ( c s )( x ) ==> ( x ) ( -s c )( y ) ( y ) where c and s are complex and the vectors x and y are complex.Parameters:NN is INTEGER The number of elements in the vectors CX and CY.CXCX is COMPLEX*16 array, dimension (N) On input, the vector x. On output, CX is overwritten with c*x + s*y.INCXINCX is INTEGER The increment between successive values of CX. INCX <> 0.CYCY is COMPLEX*16 array, dimension (N) On input, the vector y. On output, CY is overwritten with -s*x + c*y.INCYINCY is INTEGER The increment between successive values of CY. INCY <> 0.CC is COMPLEX*16SS is COMPLEX*16 C and S define the matrix [ C S ]. [ -S C ]Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016complex*16functionzladiv(complex*16X,complex*16Y)ZLADIVperforms complex division in real arithmetic, avoiding unnecessary overflow.Purpose:ZLADIV := X / Y, where X and Y are complex. The computation of X / Y will not overflow on an intermediary step unless the results overflows.Parameters:XX is COMPLEX*16YY is COMPLEX*16 The complex scalars X and Y.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaein(logicalRIGHTV,logicalNOINIT,integerN,complex*16,dimension(ldh,*)H,integerLDH,complex*16W,complex*16,dimension(*)V,complex*16,dimension(ldb,*)B,integerLDB,doubleprecision,dimension(*)RWORK,doubleprecisionEPS3,doubleprecisionSMLNUM,integerINFO)ZLAEINcomputes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.Purpose:ZLAEIN uses inverse iteration to find a right or left eigenvector corresponding to the eigenvalue W of a complex upper Hessenberg matrix H.Parameters:RIGHTVRIGHTV is LOGICAL = .TRUE. : compute right eigenvector; = .FALSE.: compute left eigenvector.NOINITNOINIT is LOGICAL = .TRUE. : no initial vector supplied in V = .FALSE.: initial vector supplied in V.NN is INTEGER The order of the matrix H. N >= 0.HH is COMPLEX*16 array, dimension (LDH,N) The upper Hessenberg matrix H.LDHLDH is INTEGER The leading dimension of the array H. LDH >= max(1,N).WW is COMPLEX*16 The eigenvalue of H whose corresponding right or left eigenvector is to be computed.VV is COMPLEX*16 array, dimension (N) On entry, if NOINIT = .FALSE., V must contain a starting vector for inverse iteration; otherwise V need not be set. On exit, V contains the computed eigenvector, normalized so that the component of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|.BB is COMPLEX*16 array, dimension (LDB,N)LDBLDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).RWORKRWORK is DOUBLE PRECISION array, dimension (N)EPS3EPS3 is DOUBLE PRECISION A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots.SMLNUMSMLNUM is DOUBLE PRECISION A machine-dependent value close to the underflow threshold.INFOINFO is INTEGER = 0: successful exit = 1: inverse iteration did not converge; V is set to the last iterate.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaev2(complex*16A,complex*16B,complex*16C,doubleprecisionRT1,doubleprecisionRT2,doubleprecisionCS1,complex*16SN1)ZLAEV2computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.Purpose:ZLAEV2 computes the eigendecomposition of a 2-by-2 Hermitian matrix [ A B ] [ CONJG(B) C ]. On return, RT1 is the eigenvalue of larger absolute value, RT2 is the eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right eigenvector for RT1, giving the decomposition [ CS1 CONJG(SN1) ] [ A B ] [ CS1 -CONJG(SN1) ] = [ RT1 0 ] [-SN1 CS1 ] [ CONJG(B) C ] [ SN1 CS1 ] [ 0 RT2 ].Parameters:AA is COMPLEX*16 The (1,1) element of the 2-by-2 matrix.BB is COMPLEX*16 The (1,2) element and the conjugate of the (2,1) element of the 2-by-2 matrix.CC is COMPLEX*16 The (2,2) element of the 2-by-2 matrix.RT1RT1 is DOUBLE PRECISION The eigenvalue of larger absolute value.RT2RT2 is DOUBLE PRECISION The eigenvalue of smaller absolute value.CS1CS1 is DOUBLE PRECISIONSN1SN1 is COMPLEX*16 The vector (CS1, SN1) is a unit right eigenvector for RT1.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:RT1 is accurate to a few ulps barring over/underflow. RT2 may be inaccurate if there is massive cancellation in the determinant A*C-B*B; higher precision or correctly rounded or correctly truncated arithmetic would be needed to compute RT2 accurately in all cases. CS1 and SN1 are accurate to a few ulps barring over/underflow. Overflow is possible only if RT1 is within a factor of 5 of overflow. Underflow is harmless if the input data is 0 or exceeds underflow_threshold / macheps.subroutinezlag2c(integerM,integerN,complex*16,dimension(lda,*)A,integerLDA,complex,dimension(ldsa,*)SA,integerLDSA,integerINFO)ZLAG2Cconverts a complex double precision matrix to a complex single precision matrix.Purpose:ZLAG2C converts a COMPLEX*16 matrix, SA, to a COMPLEX matrix, A. RMAX is the overflow for the SINGLE PRECISION arithmetic ZLAG2C checks that all the entries of A are between -RMAX and RMAX. If not the conversion is aborted and a flag is raised. This is an auxiliary routine so there is no argument checking.Parameters:MM is INTEGER The number of lines of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the M-by-N coefficient matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).SASA is COMPLEX array, dimension (LDSA,N) On exit, if INFO=0, the M-by-N coefficient matrix SA; if INFO>0, the content of SA is unspecified.LDSALDSA is INTEGER The leading dimension of the array SA. LDSA >= max(1,M).INFOINFO is INTEGER = 0: successful exit. = 1: an entry of the matrix A is greater than the SINGLE PRECISION overflow threshold, in this case, the content of SA in exit is unspecified.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlags2(logicalUPPER,doubleprecisionA1,complex*16A2,doubleprecisionA3,doubleprecisionB1,complex*16B2,doubleprecisionB3,doubleprecisionCSU,complex*16SNU,doubleprecisionCSV,complex*16SNV,doubleprecisionCSQ,complex*16SNQ)ZLAGS2Purpose:ZLAGS2 computes 2-by-2 unitary matrices U, V and Q, such that if ( UPPER ) then U**H *A*Q = U**H *( A1 A2 )*Q = ( x 0 ) ( 0 A3 ) ( x x ) and V**H*B*Q = V**H *( B1 B2 )*Q = ( x 0 ) ( 0 B3 ) ( x x ) or if ( .NOT.UPPER ) then U**H *A*Q = U**H *( A1 0 )*Q = ( x x ) ( A2 A3 ) ( 0 x ) and V**H *B*Q = V**H *( B1 0 )*Q = ( x x ) ( B2 B3 ) ( 0 x ) where U = ( CSU SNU ), V = ( CSV SNV ), ( -SNU**H CSU ) ( -SNV**H CSV ) Q = ( CSQ SNQ ) ( -SNQ**H CSQ ) The rows of the transformed A and B are parallel. Moreover, if the input 2-by-2 matrix A is not zero, then the transformed (1,1) entry of A is not zero. If the input matrices A and B are both not zero, then the transformed (2,2) element of B is not zero, except when the first rows of input A and B are parallel and the second rows are zero.Parameters:UPPERUPPER is LOGICAL = .TRUE.: the input matrices A and B are upper triangular. = .FALSE.: the input matrices A and B are lower triangular.A1A1 is DOUBLE PRECISIONA2A2 is COMPLEX*16A3A3 is DOUBLE PRECISION On entry, A1, A2 and A3 are elements of the input 2-by-2 upper (lower) triangular matrix A.B1B1 is DOUBLE PRECISIONB2B2 is COMPLEX*16B3B3 is DOUBLE PRECISION On entry, B1, B2 and B3 are elements of the input 2-by-2 upper (lower) triangular matrix B.CSUCSU is DOUBLE PRECISIONSNUSNU is COMPLEX*16 The desired unitary matrix U.CSVCSV is DOUBLE PRECISIONSNVSNV is COMPLEX*16 The desired unitary matrix V.CSQCSQ is DOUBLE PRECISIONSNQSNQ is COMPLEX*16 The desired unitary matrix Q.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlagtm(characterTRANS,integerN,integerNRHS,doubleprecisionALPHA,complex*16,dimension(*)DL,complex*16,dimension(*)D,complex*16,dimension(*)DU,complex*16,dimension(ldx,*)X,integerLDX,doubleprecisionBETA,complex*16,dimension(ldb,*)B,integerLDB)ZLAGTMperforms a matrix-matrix product of the form C = αAB+βC, where A is a tridiagonal matrix, B and C are rectangular matrices, and α and β are scalars, which may be 0, 1, or -1.Purpose:ZLAGTM performs a matrix-vector product of the form B := alpha * A * X + beta * B where A is a tridiagonal matrix of order N, B and X are N by NRHS matrices, and alpha and beta are real scalars, each of which may be 0., 1., or -1.Parameters:TRANSTRANS is CHARACTER*1 Specifies the operation applied to A. = 'N': No transpose, B := alpha * A * X + beta * B = 'T': Transpose, B := alpha * A**T * X + beta * B = 'C': Conjugate transpose, B := alpha * A**H * X + beta * BNN is INTEGER The order of the matrix A. N >= 0.NRHSNRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrices X and B.ALPHAALPHA is DOUBLE PRECISION The scalar alpha. ALPHA must be 0., 1., or -1.; otherwise, it is assumed to be 0.DLDL is COMPLEX*16 array, dimension (N-1) The (n-1) sub-diagonal elements of T.DD is COMPLEX*16 array, dimension (N) The diagonal elements of T.DUDU is COMPLEX*16 array, dimension (N-1) The (n-1) super-diagonal elements of T.XX is COMPLEX*16 array, dimension (LDX,NRHS) The N by NRHS matrix X.LDXLDX is INTEGER The leading dimension of the array X. LDX >= max(N,1).BETABETA is DOUBLE PRECISION The scalar beta. BETA must be 0., 1., or -1.; otherwise, it is assumed to be 1.BB is COMPLEX*16 array, dimension (LDB,NRHS) On entry, the N by NRHS matrix B. On exit, B is overwritten by the matrix expression B := alpha * A * X + beta * B.LDBLDB is INTEGER The leading dimension of the array B. LDB >= max(N,1).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlahqr(logicalWANTT,logicalWANTZ,integerN,integerILO,integerIHI,complex*16,dimension(ldh,*)H,integerLDH,complex*16,dimension(*)W,integerILOZ,integerIHIZ,complex*16,dimension(ldz,*)Z,integerLDZ,integerINFO)ZLAHQRcomputes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.Purpose:ZLAHQR is an auxiliary routine called by CHSEQR to update the eigenvalues and Schur decomposition already computed by CHSEQR, by dealing with the Hessenberg submatrix in rows and columns ILO to IHI.Parameters:WANTTWANTT is LOGICAL = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.WANTZWANTZ is LOGICAL = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.NN is INTEGER The order of the matrix H. N >= 0.ILOILO is INTEGERIHIIHI is INTEGER It is assumed that H is already upper triangular in rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). ZLAHQR works primarily with the Hessenberg submatrix in rows and columns ILO to IHI, but applies transformations to all of H if WANTT is .TRUE.. 1 <= ILO <= max(1,IHI); IHI <= N.HH is COMPLEX*16 array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if INFO is zero and if WANTT is .TRUE., then H is upper triangular in rows and columns ILO:IHI. If INFO is zero and if WANTT is .FALSE., then the contents of H are unspecified on exit. The output state of H in case INF is positive is below under the description of INFO.LDHLDH is INTEGER The leading dimension of the array H. LDH >= max(1,N).WW is COMPLEX*16 array, dimension (N) The computed eigenvalues ILO to IHI are stored in the corresponding elements of W. If WANTT is .TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i).ILOZILOZ is INTEGERIHIZIHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.ZZ is COMPLEX*16 array, dimension (LDZ,N) If WANTZ is .TRUE., on entry Z must contain the current matrix Z of transformations accumulated by CHSEQR, and on exit Z has been updated; transformations are applied only to the submatrix Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not referenced.LDZLDZ is INTEGER The leading dimension of the array Z. LDZ >= max(1,N).INFOINFO is INTEGER = 0: successful exit .GT. 0: if INFO = i, ZLAHQR failed to compute all the eigenvalues ILO to IHI in a total of 30 iterations per eigenvalue; elements i+1:ihi of W contain those eigenvalues which have been successfully computed. If INFO .GT. 0 and WANTT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO thorugh INFO of the final, output value of H. If INFO .GT. 0 and WANTT is .TRUE., then on exit (*) (initial value of H)*U = U*(final value of H) where U is an orthognal matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. If INFO .GT. 0 and WANTZ is .TRUE., then on exit (final value of Z) = (initial value of Z)*U where U is the orthogonal matrix in (*) (regardless of the value of WANTT.)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:02-96 Based on modifications by David Day, Sandia National Laboratory, USA 12-04 Further modifications by Ralph Byers, University of Kansas, USA This is a modified version of ZLAHQR from LAPACK version 3.0. It is (1) more robust against overflow and underflow and (2) adopts the more conservative Ahues & Tisseur stopping criterion (LAWN 122, 1997).subroutinezlahr2(integerN,integerK,integerNB,complex*16,dimension(lda,*)A,integerLDA,complex*16,dimension(nb)TAU,complex*16,dimension(ldt,nb)T,integerLDT,complex*16,dimension(ldy,nb)Y,integerLDY)ZLAHR2reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.Purpose:ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) matrix A so that elements below the k-th subdiagonal are zero. The reduction is performed by an unitary similarity transformation Q**H * A * Q. The routine returns the matrices V and T which determine Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T. This is an auxiliary routine called by ZGEHRD.Parameters:NN is INTEGER The order of the matrix A.KK is INTEGER The offset for the reduction. Elements below the k-th subdiagonal in the first NB columns are reduced to zero. K < N.NBNB is INTEGER The number of columns to be reduced.AA is COMPLEX*16 array, dimension (LDA,N-K+1) On entry, the n-by-(n-k+1) general matrix A. On exit, the elements on and above the k-th subdiagonal in the first NB columns are overwritten with the corresponding elements of the reduced matrix; the elements below the k-th subdiagonal, with the array TAU, represent the matrix Q as a product of elementary reflectors. The other columns of A are unchanged. See Further Details.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).TAUTAU is COMPLEX*16 array, dimension (NB) The scalar factors of the elementary reflectors. See Further Details.TT is COMPLEX*16 array, dimension (LDT,NB) The upper triangular matrix T.LDTLDT is INTEGER The leading dimension of the array T. LDT >= NB.YY is COMPLEX*16 array, dimension (LDY,NB) The n-by-nb matrix Y.LDYLDY is INTEGER The leading dimension of the array Y. LDY >= N.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:The matrix Q is represented as a product of nb elementary reflectors Q = H(1) H(2) . . . H(nb). Each H(i) has the form H(i) = I - tau * v * v**H where tau is a complex scalar, and v is a complex vector with v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in A(i+k+1:n,i), and tau in TAU(i). The elements of the vectors v together form the (n-k+1)-by-nb matrix V which is needed, with T and Y, to apply the transformation to the unreduced part of the matrix, using an update of the form: A := (I - V*T*V**H) * (A - Y*V**H). The contents of A on exit are illustrated by the following example with n = 7, k = 3 and nb = 2: ( a a a a a ) ( a a a a a ) ( a a a a a ) ( h h a a a ) ( v1 h a a a ) ( v1 v2 a a a ) ( v1 v2 a a a ) where a denotes an element of the original matrix A, h denotes a modified element of the upper Hessenberg matrix H, and vi denotes an element of the vector defining H(i). This subroutine is a slight modification of LAPACK-3.0's DLAHRD incorporating improvements proposed by Quintana-Orti and Van de Gejin. Note that the entries of A(1:K,2:NB) differ from those returned by the original LAPACK-3.0's DLAHRD routine. (This subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)References:Gregorio Quintana-Orti and Robert van de Geijn, 'Improving the performance of reduction to Hessenberg form,' ACM Transactions on Mathematical Software, 32(2):180-194, June 2006.subroutinezlaic1(integerJOB,integerJ,complex*16,dimension(j)X,doubleprecisionSEST,complex*16,dimension(j)W,complex*16GAMMA,doubleprecisionSESTPR,complex*16S,complex*16C)ZLAIC1applies one step of incremental condition estimation.Purpose:ZLAIC1 applies one step of incremental condition estimation in its simplest version: Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j lower triangular matrix L, such that twonorm(L*x) = sest Then ZLAIC1 computes sestpr, s, c such that the vector [ s*x ] xhat = [ c ] is an approximate singular vector of [ L 0 ] Lhat = [ w**H gamma ] in the sense that twonorm(Lhat*xhat) = sestpr. Depending on JOB, an estimate for the largest or smallest singular value is computed. Note that [s c]**H and sestpr**2 is an eigenpair of the system diag(sest*sest, 0) + [alpha gamma] * [ conjg(alpha) ] [ conjg(gamma) ] where alpha = x**H * w.Parameters:JOBJOB is INTEGER = 1: an estimate for the largest singular value is computed. = 2: an estimate for the smallest singular value is computed.JJ is INTEGER Length of X and WXX is COMPLEX*16 array, dimension (J) The j-vector x.SESTSEST is DOUBLE PRECISION Estimated singular value of j by j matrix LWW is COMPLEX*16 array, dimension (J) The j-vector w.GAMMAGAMMA is COMPLEX*16 The diagonal element gamma.SESTPRSESTPR is DOUBLE PRECISION Estimated singular value of (j+1) by (j+1) matrix Lhat.SS is COMPLEX*16 Sine needed in forming xhat.CC is COMPLEX*16 Cosine needed in forming xhat.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlangt(characterNORM,integerN,complex*16,dimension(*)DL,complex*16,dimension(*)D,complex*16,dimension(*)DU)ZLANGTreturns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general tridiagonal matrix.Purpose:ZLANGT returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex tridiagonal matrix A.Returns:ZLANGT ZLANGT = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANGT as described above.NN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANGT is set to zero.DLDL is COMPLEX*16 array, dimension (N-1) The (n-1) sub-diagonal elements of A.DD is COMPLEX*16 array, dimension (N) The diagonal elements of A.DUDU is COMPLEX*16 array, dimension (N-1) The (n-1) super-diagonal elements of A.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlanhb(characterNORM,characterUPLO,integerN,integerK,complex*16,dimension(ldab,*)AB,integerLDAB,doubleprecision,dimension(*)WORK)ZLANHBreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hermitian band matrix.Purpose:ZLANHB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n hermitian band matrix A, with k super-diagonals.Returns:ZLANHB ZLANHB = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANHB as described above.UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the band matrix A is supplied. = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANHB is set to zero.KK is INTEGER The number of super-diagonals or sub-diagonals of the band matrix A. K >= 0.ABAB is COMPLEX*16 array, dimension (LDAB,N) The upper or lower triangle of the hermitian band matrix A, stored in the first K+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+k). Note that the imaginary parts of the diagonal elements need not be set and are assumed to be zero.LDABLDAB is INTEGER The leading dimension of the array AB. LDAB >= K+1.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlanhp(characterNORM,characterUPLO,integerN,complex*16,dimension(*)AP,doubleprecision,dimension(*)WORK)ZLANHPreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form.Purpose:ZLANHP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex hermitian matrix A, supplied in packed form.Returns:ZLANHP ZLANHP = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANHP as described above.UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the hermitian matrix A is supplied. = 'U': Upper triangular part of A is supplied = 'L': Lower triangular part of A is suppliedNN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANHP is set to zero.APAP is COMPLEX*16 array, dimension (N*(N+1)/2) The upper or lower triangle of the hermitian matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. Note that the imaginary parts of the diagonal elements need not be set and are assumed to be zero.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlanhs(characterNORM,integerN,complex*16,dimension(lda,*)A,integerLDA,doubleprecision,dimension(*)WORK)ZLANHSreturns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of an upper Hessenberg matrix.Purpose:ZLANHS returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a Hessenberg matrix A.Returns:ZLANHS ZLANHS = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANHS as described above.NN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANHS is set to zero.AA is COMPLEX*16 array, dimension (LDA,N) The n by n upper Hessenberg matrix A; the part of A below the first sub-diagonal is not referenced.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(N,1).WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlanht(characterNORM,integerN,doubleprecision,dimension(*)D,complex*16,dimension(*)E)ZLANHTreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix.Purpose:ZLANHT returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian tridiagonal matrix A.Returns:ZLANHT ZLANHT = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANHT as described above.NN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANHT is set to zero.DD is DOUBLE PRECISION array, dimension (N) The diagonal elements of A.EE is COMPLEX*16 array, dimension (N-1) The (n-1) sub-diagonal or super-diagonal elements of A.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlansb(characterNORM,characterUPLO,integerN,integerK,complex*16,dimension(ldab,*)AB,integerLDAB,doubleprecision,dimension(*)WORK)ZLANSBreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix.Purpose:ZLANSB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n symmetric band matrix A, with k super-diagonals.Returns:ZLANSB ZLANSB = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANSB as described above.UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the band matrix A is supplied. = 'U': Upper triangular part is supplied = 'L': Lower triangular part is suppliedNN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANSB is set to zero.KK is INTEGER The number of super-diagonals or sub-diagonals of the band matrix A. K >= 0.ABAB is COMPLEX*16 array, dimension (LDAB,N) The upper or lower triangle of the symmetric band matrix A, stored in the first K+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+k).LDABLDAB is INTEGER The leading dimension of the array AB. LDAB >= K+1.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlansp(characterNORM,characterUPLO,integerN,complex*16,dimension(*)AP,doubleprecision,dimension(*)WORK)ZLANSPreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form.Purpose:ZLANSP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix A, supplied in packed form.Returns:ZLANSP ZLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANSP as described above.UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is supplied. = 'U': Upper triangular part of A is supplied = 'L': Lower triangular part of A is suppliedNN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANSP is set to zero.APAP is COMPLEX*16 array, dimension (N*(N+1)/2) The upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlantb(characterNORM,characterUPLO,characterDIAG,integerN,integerK,complex*16,dimension(ldab,*)AB,integerLDAB,doubleprecision,dimension(*)WORK)ZLANTBreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular band matrix.Purpose:ZLANTB returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of an n by n triangular band matrix A, with ( k + 1 ) diagonals.Returns:ZLANTB ZLANTB = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANTB as described above.UPLOUPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangularDIAGDIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangularNN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANTB is set to zero.KK is INTEGER The number of super-diagonals of the matrix A if UPLO = 'U', or the number of sub-diagonals of the matrix A if UPLO = 'L'. K >= 0.ABAB is COMPLEX*16 array, dimension (LDAB,N) The upper or lower triangular band matrix A, stored in the first k+1 rows of AB. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+k). Note that when DIAG = 'U', the elements of the array AB corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one.LDABLDAB is INTEGER The leading dimension of the array AB. LDAB >= K+1.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlantp(characterNORM,characterUPLO,characterDIAG,integerN,complex*16,dimension(*)AP,doubleprecision,dimension(*)WORK)ZLANTPreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix supplied in packed form.Purpose:ZLANTP returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a triangular matrix A, supplied in packed form.Returns:ZLANTP ZLANTP = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANTP as described above.UPLOUPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangularDIAGDIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangularNN is INTEGER The order of the matrix A. N >= 0. When N = 0, ZLANTP is set to zero.APAP is COMPLEX*16 array, dimension (N*(N+1)/2) The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. Note that when DIAG = 'U', the elements of the array AP corresponding to the diagonal elements of the matrix A are not referenced, but are assumed to be one.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= N when NORM = 'I'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016doubleprecisionfunctionzlantr(characterNORM,characterUPLO,characterDIAG,integerM,integerN,complex*16,dimension(lda,*)A,integerLDA,doubleprecision,dimension(*)WORK)ZLANTRreturns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix.Purpose:ZLANTR returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a trapezoidal or triangular matrix A.Returns:ZLANTR ZLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm' ( ( norm1(A), NORM = '1', 'O' or 'o' ( ( normI(A), NORM = 'I' or 'i' ( ( normF(A), NORM = 'F', 'f', 'E' or 'e' where norm1 denotes the one norm of a matrix (maximum column sum), normI denotes the infinity norm of a matrix (maximum row sum) and normF denotes the Frobenius norm of a matrix (square root of sum of squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.Parameters:NORMNORM is CHARACTER*1 Specifies the value to be returned in ZLANTR as described above.UPLOUPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower trapezoidal. = 'U': Upper trapezoidal = 'L': Lower trapezoidal Note that A is triangular instead of trapezoidal if M = N.DIAGDIAG is CHARACTER*1 Specifies whether or not the matrix A has unit diagonal. = 'N': Non-unit diagonal = 'U': Unit diagonalMM is INTEGER The number of rows of the matrix A. M >= 0, and if UPLO = 'U', M <= N. When M = 0, ZLANTR is set to zero.NN is INTEGER The number of columns of the matrix A. N >= 0, and if UPLO = 'L', N <= M. When N = 0, ZLANTR is set to zero.AA is COMPLEX*16 array, dimension (LDA,N) The trapezoidal matrix A (A is triangular if M = N). If UPLO = 'U', the leading m by n upper trapezoidal part of the array A contains the upper trapezoidal matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading m by n lower trapezoidal part of the array A contains the lower trapezoidal matrix, and the strictly upper triangular part of A is not referenced. Note that when DIAG = 'U', the diagonal elements of A are not referenced and are assumed to be one.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(M,1).WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)), where LWORK >= M when NORM = 'I'; otherwise, WORK is not referenced.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlapll(integerN,complex*16,dimension(*)X,integerINCX,complex*16,dimension(*)Y,integerINCY,doubleprecisionSSMIN)ZLAPLLmeasures the linear dependence of two vectors.Purpose:Given two column vectors X and Y, let A = ( X Y ). The subroutine first computes the QR factorization of A = Q*R, and then computes the SVD of the 2-by-2 upper triangular matrix R. The smaller singular value of R is returned in SSMIN, which is used as the measurement of the linear dependency of the vectors X and Y.Parameters:NN is INTEGER The length of the vectors X and Y.XX is COMPLEX*16 array, dimension (1+(N-1)*INCX) On entry, X contains the N-vector X. On exit, X is overwritten.INCXINCX is INTEGER The increment between successive elements of X. INCX > 0.YY is COMPLEX*16 array, dimension (1+(N-1)*INCY) On entry, Y contains the N-vector Y. On exit, Y is overwritten.INCYINCY is INTEGER The increment between successive elements of Y. INCY > 0.SSMINSSMIN is DOUBLE PRECISION The smallest singular value of the N-by-2 matrix A = ( X Y ).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlapmr(logicalFORWRD,integerM,integerN,complex*16,dimension(ldx,*)X,integerLDX,integer,dimension(*)K)ZLAPMRrearranges rows of a matrix as specified by a permutation vector.Purpose:ZLAPMR rearranges the rows of the M by N matrix X as specified by the permutation K(1),K(2),...,K(M) of the integers 1,...,M. If FORWRD = .TRUE., forward permutation: X(K(I),*) is moved X(I,*) for I = 1,2,...,M. If FORWRD = .FALSE., backward permutation: X(I,*) is moved to X(K(I),*) for I = 1,2,...,M.Parameters:FORWRDFORWRD is LOGICAL = .TRUE., forward permutation = .FALSE., backward permutationMM is INTEGER The number of rows of the matrix X. M >= 0.NN is INTEGER The number of columns of the matrix X. N >= 0.XX is COMPLEX*16 array, dimension (LDX,N) On entry, the M by N matrix X. On exit, X contains the permuted matrix X.LDXLDX is INTEGER The leading dimension of the array X, LDX >= MAX(1,M).KK is INTEGER array, dimension (M) On entry, K contains the permutation vector. K is used as internal workspace, but reset to its original value on output.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlapmt(logicalFORWRD,integerM,integerN,complex*16,dimension(ldx,*)X,integerLDX,integer,dimension(*)K)ZLAPMTperforms a forward or backward permutation of the columns of a matrix.Purpose:ZLAPMT rearranges the columns of the M by N matrix X as specified by the permutation K(1),K(2),...,K(N) of the integers 1,...,N. If FORWRD = .TRUE., forward permutation: X(*,K(J)) is moved X(*,J) for J = 1,2,...,N. If FORWRD = .FALSE., backward permutation: X(*,J) is moved to X(*,K(J)) for J = 1,2,...,N.Parameters:FORWRDFORWRD is LOGICAL = .TRUE., forward permutation = .FALSE., backward permutationMM is INTEGER The number of rows of the matrix X. M >= 0.NN is INTEGER The number of columns of the matrix X. N >= 0.XX is COMPLEX*16 array, dimension (LDX,N) On entry, the M by N matrix X. On exit, X contains the permuted matrix X.LDXLDX is INTEGER The leading dimension of the array X, LDX >= MAX(1,M).KK is INTEGER array, dimension (N) On entry, K contains the permutation vector. K is used as internal workspace, but reset to its original value on output.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaqhb(characterUPLO,integerN,integerKD,complex*16,dimension(ldab,*)AB,integerLDAB,doubleprecision,dimension(*)S,doubleprecisionSCOND,doubleprecisionAMAX,characterEQUED)ZLAQHBscales a Hermitian band matrix, using scaling factors computed by cpbequ.Purpose:ZLAQHB equilibrates a Hermitian band matrix A using the scaling factors in the vector S.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored. = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the matrix A. N >= 0.KDKD 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.ABAB is COMPLEX*16 array, dimension (LDAB,N) On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**H *U or A = L*L**H of the band matrix A, in the same storage format as A.LDABLDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.SS is DOUBLE PRECISION array, dimension (N) The scale factors for A.SCONDSCOND is DOUBLE PRECISION Ratio of the smallest S(i) to the largest S(i).AMAXAMAX is DOUBLE PRECISION Absolute value of largest matrix entry.EQUEDEQUED is CHARACTER*1 Specifies whether or not equilibration was done. = 'N': No equilibration. = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).InternalParameters:THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done. LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaqhp(characterUPLO,integerN,complex*16,dimension(*)AP,doubleprecision,dimension(*)S,doubleprecisionSCOND,doubleprecisionAMAX,characterEQUED)ZLAQHPscales a Hermitian matrix stored in packed form.Purpose:ZLAQHP equilibrates a Hermitian matrix A using the scaling factors in the vector S.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the Hermitian matrix A is stored. = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the matrix A. N >= 0.APAP is COMPLEX*16 array, dimension (N*(N+1)/2) On entry, the upper or lower triangle of the Hermitian matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. On exit, the equilibrated matrix: diag(S) * A * diag(S), in the same storage format as A.SS is DOUBLE PRECISION array, dimension (N) The scale factors for A.SCONDSCOND is DOUBLE PRECISION Ratio of the smallest S(i) to the largest S(i).AMAXAMAX is DOUBLE PRECISION Absolute value of largest matrix entry.EQUEDEQUED is CHARACTER*1 Specifies whether or not equilibration was done. = 'N': No equilibration. = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).InternalParameters:THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done. LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaqp2(integerM,integerN,integerOFFSET,complex*16,dimension(lda,*)A,integerLDA,integer,dimension(*)JPVT,complex*16,dimension(*)TAU,doubleprecision,dimension(*)VN1,doubleprecision,dimension(*)VN2,complex*16,dimension(*)WORK)ZLAQP2computes a QR factorization with column pivoting of the matrix block.Purpose:ZLAQP2 computes a QR factorization with column pivoting of the block A(OFFSET+1:M,1:N). The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.Parameters:MM is INTEGER The number of rows of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0.OFFSETOFFSET is INTEGER The number of rows of the matrix A that must be pivoted but no factorized. OFFSET >= 0.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, the upper triangle of block A(OFFSET+1:M,1:N) is the triangular factor obtained; the elements in block A(OFFSET+1:M,1:N) below the diagonal, together with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors. Block A(1:OFFSET,1:N) has been accordingly pivoted, but no factorized.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).JPVTJPVT is INTEGER array, dimension (N) On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted to the front of A*P (a leading column); if JPVT(i) = 0, the i-th column of A is a free column. On exit, if JPVT(i) = k, then the i-th column of A*P was the k-th column of A.TAUTAU is COMPLEX*16 array, dimension (min(M,N)) The scalar factors of the elementary reflectors.VN1VN1 is DOUBLE PRECISION array, dimension (N) The vector with the partial column norms.VN2VN2 is DOUBLE PRECISION array, dimension (N) The vector with the exact column norms.WORKWORK is COMPLEX*16 array, dimension (N)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.References:LAPACK Working Note 176subroutinezlaqps(integerM,integerN,integerOFFSET,integerNB,integerKB,complex*16,dimension(lda,*)A,integerLDA,integer,dimension(*)JPVT,complex*16,dimension(*)TAU,doubleprecision,dimension(*)VN1,doubleprecision,dimension(*)VN2,complex*16,dimension(*)AUXV,complex*16,dimension(ldf,*)F,integerLDF)ZLAQPScomputes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3.Purpose:ZLAQPS computes a step of QR factorization with column pivoting of a complex M-by-N matrix A by using Blas-3. It tries to factorize NB columns from A starting from the row OFFSET+1, and updates all of the matrix with Blas-3 xGEMM. In some cases, due to catastrophic cancellations, it cannot factorize NB columns. Hence, the actual number of factorized columns is returned in KB. Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.Parameters:MM is INTEGER The number of rows of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0OFFSETOFFSET is INTEGER The number of rows of A that have been factorized in previous steps.NBNB is INTEGER The number of columns to factorize.KBKB is INTEGER The number of columns actually factorized.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, block A(OFFSET+1:M,1:KB) is the triangular factor obtained and block A(1:OFFSET,1:N) has been accordingly pivoted, but no factorized. The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has been updated.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).JPVTJPVT is INTEGER array, dimension (N) JPVT(I) = K <==> Column K of the full matrix A has been permuted into position I in AP.TAUTAU is COMPLEX*16 array, dimension (KB) The scalar factors of the elementary reflectors.VN1VN1 is DOUBLE PRECISION array, dimension (N) The vector with the partial column norms.VN2VN2 is DOUBLE PRECISION array, dimension (N) The vector with the exact column norms.AUXVAUXV is COMPLEX*16 array, dimension (NB) Auxiliar vector.FF is COMPLEX*16 array, dimension (LDF,NB) Matrix F**H = L * Y**H * A.LDFLDF is INTEGER The leading dimension of the array F. LDF >= max(1,N).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.References:LAPACK Working Note 176subroutinezlaqr0(logicalWANTT,logicalWANTZ,integerN,integerILO,integerIHI,complex*16,dimension(ldh,*)H,integerLDH,complex*16,dimension(*)W,integerILOZ,integerIHIZ,complex*16,dimension(ldz,*)Z,integerLDZ,complex*16,dimension(*)WORK,integerLWORK,integerINFO)ZLAQR0computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.Purpose:ZLAQR0 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**H, where T is an upper triangular matrix (the Schur form), and Z is the unitary matrix of Schur vectors. Optionally Z may be postmultiplied into an input unitary matrix Q so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.Parameters:WANTTWANTT is LOGICAL = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.WANTZWANTZ is LOGICAL = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.NN is INTEGER The order of the matrix H. N .GE. 0.ILOILO is INTEGERIHIIHI is INTEGER It is assumed that H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, H(ILO,ILO-1) is zero. ILO and IHI are normally set by a previous call to ZGEBAL, and then passed to ZGEHRD when the matrix output by ZGEBAL is reduced to Hessenberg form. Otherwise, ILO and IHI should be set to 1 and N, respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. If N = 0, then ILO = 1 and IHI = 0.HH is COMPLEX*16 array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper triangular matrix T from the Schur decomposition (the Schur form). If INFO = 0 and WANT is .FALSE., then the contents of H are unspecified on exit. (The output value of H when INFO.GT.0 is given under the description of INFO below.) This subroutine may explicitly set H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.LDHLDH is INTEGER The leading dimension of the array H. LDH .GE. max(1,N).WW is COMPLEX*16 array, dimension (N) The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i).ILOZILOZ is INTEGERIHIZIHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.ZZ is COMPLEX*16 array, dimension (LDZ,IHI) If WANTZ is .FALSE., then Z is not referenced. If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the orthogonal Schur factor of H(ILO:IHI,ILO:IHI). (The output value of Z when INFO.GT.0 is given under the description of INFO below.)LDZLDZ is INTEGER The leading dimension of the array Z. if WANTZ is .TRUE. then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.WORKWORK is COMPLEX*16 array, dimension LWORK On exit, if LWORK = -1, WORK(1) returns an estimate of the optimal value for LWORK.LWORKLWORK is INTEGER The dimension of the array WORK. LWORK .GE. max(1,N) is sufficient, but LWORK typically as large as 6*N may be required for optimal performance. A workspace query to determine the optimal workspace size is recommended. If LWORK = -1, then ZLAQR0 does a workspace query. In this case, ZLAQR0 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.INFOINFO is INTEGER = 0: successful exit .GT. 0: if INFO = i, ZLAQR0 failed to compute all of the eigenvalues. Elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.) If INFO .GT. 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigen- values of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. If INFO .GT. 0 and WANTT is .TRUE., then on exit (*) (initial value of H)*U = U*(final value of H) where U is a unitary matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. If INFO .GT. 0 and WANTZ is .TRUE., then on exit (final value of Z(ILO:IHI,ILOZ:IHIZ) = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U where U is the unitary matrix in (*) (regard- less of the value of WANTT.) If INFO .GT. 0 and WANTZ is .FALSE., then Z is not accessed.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USAReferences: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.subroutinezlaqr1(integerN,complex*16,dimension(ldh,*)H,integerLDH,complex*16S1,complex*16S2,complex*16,dimension(*)V)ZLAQR1sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and specified shifts.Purpose:Given a 2-by-2 or 3-by-3 matrix H, ZLAQR1 sets v to a scalar multiple of the first column of the product (*) K = (H - s1*I)*(H - s2*I) scaling to avoid overflows and most underflows. This is useful for starting double implicit shift bulges in the QR algorithm.Parameters:NN is INTEGER Order of the matrix H. N must be either 2 or 3.HH is COMPLEX*16 array, dimension (LDH,N) The 2-by-2 or 3-by-3 matrix H in (*).LDHLDH is INTEGER The leading dimension of H as declared in the calling procedure. LDH.GE.NS1S1 is COMPLEX*16S2S2 is COMPLEX*16 S1 and S2 are the shifts defining K in (*) above.VV is COMPLEX*16 array, dimension (N) A scalar multiple of the first column of the matrix K in (*).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017Contributors:Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USAsubroutinezlaqr2(logicalWANTT,logicalWANTZ,integerN,integerKTOP,integerKBOT,integerNW,complex*16,dimension(ldh,*)H,integerLDH,integerILOZ,integerIHIZ,complex*16,dimension(ldz,*)Z,integerLDZ,integerNS,integerND,complex*16,dimension(*)SH,complex*16,dimension(ldv,*)V,integerLDV,integerNH,complex*16,dimension(ldt,*)T,integerLDT,integerNV,complex*16,dimension(ldwv,*)WV,integerLDWV,complex*16,dimension(*)WORK,integerLWORK)ZLAQR2performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).Purpose:ZLAQR2 is identical to ZLAQR3 except that it avoids recursion by calling ZLAHQR instead of ZLAQR4. Aggressive early deflation: ZLAQR2 accepts as input an upper Hessenberg matrix H and performs an unitary similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been over- written by a new Hessenberg matrix that is a perturbation of an unitary similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries.Parameters:WANTTWANTT is LOGICAL If .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues.WANTZWANTZ is LOGICAL If .TRUE., then the unitary matrix Z is updated so so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced.NN is INTEGER The order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary matrix Z.KTOPKTOP is INTEGER It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.KBOTKBOT is INTEGER It is assumed without a check that either KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.NWNW is INTEGER Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).HH is COMPLEX*16 array, dimension (LDH,N) On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries.LDHLDH is INTEGER Leading dimension of H just as declared in the calling subroutine. N .LE. LDHILOZILOZ is INTEGERIHIZIHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.ZZ is COMPLEX*16 array, dimension (LDZ,N) IF WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ is .FALSE., then Z is unreferenced.LDZLDZ is INTEGER The leading dimension of Z just as declared in the calling subroutine. 1 .LE. LDZ.NSNS is INTEGER The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine.NDND is INTEGER The number of converged eigenvalues uncovered by this subroutine.SHSH is COMPLEX*16 array, dimension (KBOT) On output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT-ND). Converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT).VV is COMPLEX*16 array, dimension (LDV,NW) An NW-by-NW work array.LDVLDV is INTEGER The leading dimension of V just as declared in the calling subroutine. NW .LE. LDVNHNH is INTEGER The number of columns of T. NH.GE.NW.TT is COMPLEX*16 array, dimension (LDT,NW)LDTLDT is INTEGER The leading dimension of T just as declared in the calling subroutine. NW .LE. LDTNVNV is INTEGER The number of rows of work array WV available for workspace. NV.GE.NW.WVWV is COMPLEX*16 array, dimension (LDWV,NW)LDWVLDWV is INTEGER The leading dimension of W just as declared in the calling subroutine. NW .LE. LDVWORKWORK is COMPLEX*16 array, dimension (LWORK) On exit, WORK(1) is set to an estimate of the optimal value of LWORK for the given values of N, NW, KTOP and KBOT.LWORKLWORK is INTEGER The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK. If LWORK = -1, then a workspace query is assumed; ZLAQR2 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017Contributors:Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USAsubroutinezlaqr3(logicalWANTT,logicalWANTZ,integerN,integerKTOP,integerKBOT,integerNW,complex*16,dimension(ldh,*)H,integerLDH,integerILOZ,integerIHIZ,complex*16,dimension(ldz,*)Z,integerLDZ,integerNS,integerND,complex*16,dimension(*)SH,complex*16,dimension(ldv,*)V,integerLDV,integerNH,complex*16,dimension(ldt,*)T,integerLDT,integerNV,complex*16,dimension(ldwv,*)WV,integerLDWV,complex*16,dimension(*)WORK,integerLWORK)ZLAQR3performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).Purpose:Aggressive early deflation: ZLAQR3 accepts as input an upper Hessenberg matrix H and performs an unitary similarity transformation designed to detect and deflate fully converged eigenvalues from a trailing principal submatrix. On output H has been over- written by a new Hessenberg matrix that is a perturbation of an unitary similarity transformation of H. It is to be hoped that the final version of H has many zero subdiagonal entries.Parameters:WANTTWANTT is LOGICAL If .TRUE., then the Hessenberg matrix H is fully updated so that the triangular Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then only enough of H is updated to preserve the eigenvalues.WANTZWANTZ is LOGICAL If .TRUE., then the unitary matrix Z is updated so so that the unitary Schur factor may be computed (in cooperation with the calling subroutine). If .FALSE., then Z is not referenced.NN is INTEGER The order of the matrix H and (if WANTZ is .TRUE.) the order of the unitary matrix Z.KTOPKTOP is INTEGER It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.KBOTKBOT is INTEGER It is assumed without a check that either KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together determine an isolated block along the diagonal of the Hessenberg matrix.NWNW is INTEGER Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).HH is COMPLEX*16 array, dimension (LDH,N) On input the initial N-by-N section of H stores the Hessenberg matrix undergoing aggressive early deflation. On output H has been transformed by a unitary similarity transformation, perturbed, and the returned to Hessenberg form that (it is to be hoped) has some zero subdiagonal entries.LDHLDH is INTEGER Leading dimension of H just as declared in the calling subroutine. N .LE. LDHILOZILOZ is INTEGERIHIZIHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.ZZ is COMPLEX*16 array, dimension (LDZ,N) IF WANTZ is .TRUE., then on output, the unitary similarity transformation mentioned above has been accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ is .FALSE., then Z is unreferenced.LDZLDZ is INTEGER The leading dimension of Z just as declared in the calling subroutine. 1 .LE. LDZ.NSNS is INTEGER The number of unconverged (ie approximate) eigenvalues returned in SR and SI that may be used as shifts by the calling subroutine.NDND is INTEGER The number of converged eigenvalues uncovered by this subroutine.SHSH is COMPLEX*16 array, dimension (KBOT) On output, approximate eigenvalues that may be used for shifts are stored in SH(KBOT-ND-NS+1) through SR(KBOT-ND). Converged eigenvalues are stored in SH(KBOT-ND+1) through SH(KBOT).VV is COMPLEX*16 array, dimension (LDV,NW) An NW-by-NW work array.LDVLDV is INTEGER The leading dimension of V just as declared in the calling subroutine. NW .LE. LDVNHNH is INTEGER The number of columns of T. NH.GE.NW.TT is COMPLEX*16 array, dimension (LDT,NW)LDTLDT is INTEGER The leading dimension of T just as declared in the calling subroutine. NW .LE. LDTNVNV is INTEGER The number of rows of work array WV available for workspace. NV.GE.NW.WVWV is COMPLEX*16 array, dimension (LDWV,NW)LDWVLDWV is INTEGER The leading dimension of W just as declared in the calling subroutine. NW .LE. LDVWORKWORK is COMPLEX*16 array, dimension (LWORK) On exit, WORK(1) is set to an estimate of the optimal value of LWORK for the given values of N, NW, KTOP and KBOT.LWORKLWORK is INTEGER The dimension of the work array WORK. LWORK = 2*NW suffices, but greater efficiency may result from larger values of LWORK. If LWORK = -1, then a workspace query is assumed; ZLAQR3 only estimates the optimal workspace size for the given values of N, NW, KTOP and KBOT. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USAsubroutinezlaqr4(logicalWANTT,logicalWANTZ,integerN,integerILO,integerIHI,complex*16,dimension(ldh,*)H,integerLDH,complex*16,dimension(*)W,integerILOZ,integerIHIZ,complex*16,dimension(ldz,*)Z,integerLDZ,complex*16,dimension(*)WORK,integerLWORK,integerINFO)ZLAQR4computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.Purpose:ZLAQR4 implements one level of recursion for ZLAQR0. It is a complete implementation of the small bulge multi-shift QR algorithm. It may be called by ZLAQR0 and, for large enough deflation window size, it may be called by ZLAQR3. This subroutine is identical to ZLAQR0 except that it calls ZLAQR2 instead of ZLAQR3. ZLAQR4 computes the eigenvalues of a Hessenberg matrix H and, optionally, the matrices T and Z from the Schur decomposition H = Z T Z**H, where T is an upper triangular matrix (the Schur form), and Z is the unitary matrix of Schur vectors. Optionally Z may be postmultiplied into an input unitary matrix Q so that this routine can give the Schur factorization of a matrix A which has been reduced to the Hessenberg form H by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.Parameters:WANTTWANTT is LOGICAL = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.WANTZWANTZ is LOGICAL = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.NN is INTEGER The order of the matrix H. N .GE. 0.ILOILO is INTEGERIHIIHI is INTEGER It is assumed that H is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, H(ILO,ILO-1) is zero. ILO and IHI are normally set by a previous call to ZGEBAL, and then passed to ZGEHRD when the matrix output by ZGEBAL is reduced to Hessenberg form. Otherwise, ILO and IHI should be set to 1 and N, respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. If N = 0, then ILO = 1 and IHI = 0.HH is COMPLEX*16 array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if INFO = 0 and WANTT is .TRUE., then H contains the upper triangular matrix T from the Schur decomposition (the Schur form). If INFO = 0 and WANT is .FALSE., then the contents of H are unspecified on exit. (The output value of H when INFO.GT.0 is given under the description of INFO below.) This subroutine may explicitly set H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.LDHLDH is INTEGER The leading dimension of the array H. LDH .GE. max(1,N).WW is COMPLEX*16 array, dimension (N) The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with W(i) = H(i,i).ILOZILOZ is INTEGERIHIZIHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.ZZ is COMPLEX*16 array, dimension (LDZ,IHI) If WANTZ is .FALSE., then Z is not referenced. If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the orthogonal Schur factor of H(ILO:IHI,ILO:IHI). (The output value of Z when INFO.GT.0 is given under the description of INFO below.)LDZLDZ is INTEGER The leading dimension of the array Z. if WANTZ is .TRUE. then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.WORKWORK is COMPLEX*16 array, dimension LWORK On exit, if LWORK = -1, WORK(1) returns an estimate of the optimal value for LWORK.LWORKLWORK is INTEGER The dimension of the array WORK. LWORK .GE. max(1,N) is sufficient, but LWORK typically as large as 6*N may be required for optimal performance. A workspace query to determine the optimal workspace size is recommended. If LWORK = -1, then ZLAQR4 does a workspace query. In this case, ZLAQR4 checks the input parameters and estimates the optimal workspace size for the given values of N, ILO and IHI. The estimate is returned in WORK(1). No error message related to LWORK is issued by XERBLA. Neither H nor Z are accessed.INFOINFO is INTEGER = 0: successful exit .GT. 0: if INFO = i, ZLAQR4 failed to compute all of the eigenvalues. Elements 1:ilo-1 and i+1:n of WR and WI contain those eigenvalues which have been successfully computed. (Failures are rare.) If INFO .GT. 0 and WANT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigen- values of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. If INFO .GT. 0 and WANTT is .TRUE., then on exit (*) (initial value of H)*U = U*(final value of H) where U is a unitary matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. If INFO .GT. 0 and WANTZ is .TRUE., then on exit (final value of Z(ILO:IHI,ILOZ:IHIZ) = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U where U is the unitary matrix in (*) (regard- less of the value of WANTT.) If INFO .GT. 0 and WANTZ is .FALSE., then Z is not accessed.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USAReferences: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.subroutinezlaqr5(logicalWANTT,logicalWANTZ,integerKACC22,integerN,integerKTOP,integerKBOT,integerNSHFTS,complex*16,dimension(*)S,complex*16,dimension(ldh,*)H,integerLDH,integerILOZ,integerIHIZ,complex*16,dimension(ldz,*)Z,integerLDZ,complex*16,dimension(ldv,*)V,integerLDV,complex*16,dimension(ldu,*)U,integerLDU,integerNV,complex*16,dimension(ldwv,*)WV,integerLDWV,integerNH,complex*16,dimension(ldwh,*)WH,integerLDWH)ZLAQR5performs a single small-bulge multi-shift QR sweep.Purpose:ZLAQR5, called by ZLAQR0, performs a single small-bulge multi-shift QR sweep.Parameters:WANTTWANTT is LOGICAL WANTT = .true. if the triangular Schur factor is being computed. WANTT is set to .false. otherwise.WANTZWANTZ is LOGICAL WANTZ = .true. if the unitary Schur factor is being computed. WANTZ is set to .false. otherwise.KACC22KACC22 is INTEGER with value 0, 1, or 2. Specifies the computation mode of far-from-diagonal orthogonal updates. = 0: ZLAQR5 does not accumulate reflections and does not use matrix-matrix multiply to update far-from-diagonal matrix entries. = 1: ZLAQR5 accumulates reflections and uses matrix-matrix multiply to update the far-from-diagonal matrix entries. = 2: ZLAQR5 accumulates reflections, uses matrix-matrix multiply to update the far-from-diagonal matrix entries, and takes advantage of 2-by-2 block structure during matrix multiplies.NN is INTEGER N is the order of the Hessenberg matrix H upon which this subroutine operates.KTOPKTOP is INTEGERKBOTKBOT is INTEGER These are the first and last rows and columns of an isolated diagonal block upon which the QR sweep is to be applied. It is assumed without a check that either KTOP = 1 or H(KTOP,KTOP-1) = 0 and either KBOT = N or H(KBOT+1,KBOT) = 0.NSHFTSNSHFTS is INTEGER NSHFTS gives the number of simultaneous shifts. NSHFTS must be positive and even.SS is COMPLEX*16 array, dimension (NSHFTS) S contains the shifts of origin that define the multi- shift QR sweep. On output S may be reordered.HH is COMPLEX*16 array, dimension (LDH,N) On input H contains a Hessenberg matrix. On output a multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied to the isolated diagonal block in rows and columns KTOP through KBOT.LDHLDH is INTEGER LDH is the leading dimension of H just as declared in the calling procedure. LDH.GE.MAX(1,N).ILOZILOZ is INTEGERIHIZIHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. NZZ is COMPLEX*16 array, dimension (LDZ,IHIZ) If WANTZ = .TRUE., then the QR Sweep unitary similarity transformation is accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. If WANTZ = .FALSE., then Z is unreferenced.LDZLDZ is INTEGER LDA is the leading dimension of Z just as declared in the calling procedure. LDZ.GE.N.VV is COMPLEX*16 array, dimension (LDV,NSHFTS/2)LDVLDV is INTEGER LDV is the leading dimension of V as declared in the calling procedure. LDV.GE.3.UU is COMPLEX*16 array, dimension (LDU,3*NSHFTS-3)LDULDU is INTEGER LDU is the leading dimension of U just as declared in the in the calling subroutine. LDU.GE.3*NSHFTS-3.NHNH is INTEGER NH is the number of columns in array WH available for workspace. NH.GE.1.WHWH is COMPLEX*16 array, dimension (LDWH,NH)LDWHLDWH is INTEGER Leading dimension of WH just as declared in the calling procedure. LDWH.GE.3*NSHFTS-3.NVNV is INTEGER NV is the number of rows in WV agailable for workspace. NV.GE.1.WVWV is COMPLEX*16 array, dimension (LDWV,3*NSHFTS-3)LDWVLDWV is INTEGER LDWV is the leading dimension of WV as declared in the in the calling subroutine. LDWV.GE.NV.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USAReferences: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.subroutinezlaqsb(characterUPLO,integerN,integerKD,complex*16,dimension(ldab,*)AB,integerLDAB,doubleprecision,dimension(*)S,doubleprecisionSCOND,doubleprecisionAMAX,characterEQUED)ZLAQSBscales a symmetric/Hermitian band matrix, using scaling factors computed by spbequ.Purpose:ZLAQSB equilibrates a symmetric band matrix A using the scaling factors in the vector S.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored. = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the matrix A. N >= 0.KDKD 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.ABAB is COMPLEX*16 array, dimension (LDAB,N) On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**H *U or A = L*L**H of the band matrix A, in the same storage format as A.LDABLDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.SS is DOUBLE PRECISION array, dimension (N) The scale factors for A.SCONDSCOND is DOUBLE PRECISION Ratio of the smallest S(i) to the largest S(i).AMAXAMAX is DOUBLE PRECISION Absolute value of largest matrix entry.EQUEDEQUED is CHARACTER*1 Specifies whether or not equilibration was done. = 'N': No equilibration. = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).InternalParameters:THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done. LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaqsp(characterUPLO,integerN,complex*16,dimension(*)AP,doubleprecision,dimension(*)S,doubleprecisionSCOND,doubleprecisionAMAX,characterEQUED)ZLAQSPscales a symmetric/Hermitian matrix in packed storage, using scaling factors computed by sppequ.Purpose:ZLAQSP equilibrates a symmetric matrix A using the scaling factors in the vector S.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored. = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the matrix A. N >= 0.APAP is COMPLEX*16 array, dimension (N*(N+1)/2) On entry, the upper or lower triangle of the symmetric matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. On exit, the equilibrated matrix: diag(S) * A * diag(S), in the same storage format as A.SS is DOUBLE PRECISION array, dimension (N) The scale factors for A.SCONDSCOND is DOUBLE PRECISION Ratio of the smallest S(i) to the largest S(i).AMAXAMAX is DOUBLE PRECISION Absolute value of largest matrix entry.EQUEDEQUED is CHARACTER*1 Specifies whether or not equilibration was done. = 'N': No equilibration. = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S).InternalParameters:THRESH is a threshold value used to decide if scaling should be done based on the ratio of the scaling factors. If SCOND < THRESH, scaling is done. LARGE and SMALL are threshold values used to decide if scaling should be done based on the absolute size of the largest matrix element. If AMAX > LARGE or AMAX < SMALL, scaling is done.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlar1v(integerN,integerB1,integerBN,doubleprecisionLAMBDA,doubleprecision,dimension(*)D,doubleprecision,dimension(*)L,doubleprecision,dimension(*)LD,doubleprecision,dimension(*)LLD,doubleprecisionPIVMIN,doubleprecisionGAPTOL,complex*16,dimension(*)Z,logicalWANTNC,integerNEGCNT,doubleprecisionZTZ,doubleprecisionMINGMA,integerR,integer,dimension(*)ISUPPZ,doubleprecisionNRMINV,doubleprecisionRESID,doubleprecisionRQCORR,doubleprecision,dimension(*)WORK)ZLAR1Vcomputes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.Purpose:ZLAR1V computes the (scaled) r-th column of the inverse of the sumbmatrix in rows B1 through BN of the tridiagonal matrix L D L**T - sigma I. When sigma is close to an eigenvalue, the computed vector is an accurate eigenvector. Usually, r corresponds to the index where the eigenvector is largest in magnitude. The following steps accomplish this computation : (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T, (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T, (c) Computation of the diagonal elements of the inverse of L D L**T - sigma I by combining the above transforms, and choosing r as the index where the diagonal of the inverse is (one of the) largest in magnitude. (d) Computation of the (scaled) r-th column of the inverse using the twisted factorization obtained by combining the top part of the the stationary and the bottom part of the progressive transform.Parameters:NN is INTEGER The order of the matrix L D L**T.B1B1 is INTEGER First index of the submatrix of L D L**T.BNBN is INTEGER Last index of the submatrix of L D L**T.LAMBDALAMBDA is DOUBLE PRECISION The shift. In order to compute an accurate eigenvector, LAMBDA should be a good approximation to an eigenvalue of L D L**T.LL is DOUBLE PRECISION array, dimension (N-1) The (n-1) subdiagonal elements of the unit bidiagonal matrix L, in elements 1 to N-1.DD is DOUBLE PRECISION array, dimension (N) The n diagonal elements of the diagonal matrix D.LDLD is DOUBLE PRECISION array, dimension (N-1) The n-1 elements L(i)*D(i).LLDLLD is DOUBLE PRECISION array, dimension (N-1) The n-1 elements L(i)*L(i)*D(i).PIVMINPIVMIN is DOUBLE PRECISION The minimum pivot in the Sturm sequence.GAPTOLGAPTOL is DOUBLE PRECISION Tolerance that indicates when eigenvector entries are negligible w.r.t. their contribution to the residual.ZZ is COMPLEX*16 array, dimension (N) On input, all entries of Z must be set to 0. On output, Z contains the (scaled) r-th column of the inverse. The scaling is such that Z(R) equals 1.WANTNCWANTNC is LOGICAL Specifies whether NEGCNT has to be computed.NEGCNTNEGCNT is INTEGER If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.ZTZZTZ is DOUBLE PRECISION The square of the 2-norm of Z.MINGMAMINGMA is DOUBLE PRECISION The reciprocal of the largest (in magnitude) diagonal element of the inverse of L D L**T - sigma I.RR is INTEGER The twist index for the twisted factorization used to compute Z. On input, 0 <= R <= N. If R is input as 0, R is set to the index where (L D L**T - sigma I)^{-1} is largest in magnitude. If 1 <= R <= N, R is unchanged. On output, R contains the twist index used to compute Z. Ideally, R designates the position of the maximum entry in the eigenvector.ISUPPZISUPPZ is INTEGER array, dimension (2) The support of the vector in Z, i.e., the vector Z is nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).NRMINVNRMINV is DOUBLE PRECISION NRMINV = 1/SQRT( ZTZ )RESIDRESID is DOUBLE PRECISION The residual of the FP vector. RESID = ABS( MINGMA )/SQRT( ZTZ )RQCORRRQCORR is DOUBLE PRECISION The Rayleigh Quotient correction to LAMBDA. RQCORR = MINGMA*TMPWORKWORK is DOUBLE PRECISION array, dimension (4*N)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Beresford Parlett, University of California, Berkeley, USA Jim Demmel, University of California, Berkeley, USA Inderjit Dhillon, University of Texas, Austin, USA Osni Marques, LBNL/NERSC, USA Christof Voemel, University of California, Berkeley, USAsubroutinezlar2v(integerN,complex*16,dimension(*)X,complex*16,dimension(*)Y,complex*16,dimension(*)Z,integerINCX,doubleprecision,dimension(*)C,complex*16,dimension(*)S,integerINCC)ZLAR2Vapplies a vector of plane rotations with real cosines and complex sines from both sides to a sequence of 2-by-2 symmetric/Hermitian matrices.Purpose:ZLAR2V applies a vector of complex plane rotations with real cosines from both sides to a sequence of 2-by-2 complex Hermitian matrices, defined by the elements of the vectors x, y and z. For i = 1,2,...,n ( x(i) z(i) ) := ( conjg(z(i)) y(i) ) ( c(i) conjg(s(i)) ) ( x(i) z(i) ) ( c(i) -conjg(s(i)) ) ( -s(i) c(i) ) ( conjg(z(i)) y(i) ) ( s(i) c(i) )Parameters:NN is INTEGER The number of plane rotations to be applied.XX is COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector x; the elements of x are assumed to be real.YY is COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector y; the elements of y are assumed to be real.ZZ is COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector z.INCXINCX is INTEGER The increment between elements of X, Y and Z. INCX > 0.CC is DOUBLE PRECISION array, dimension (1+(N-1)*INCC) The cosines of the plane rotations.SS is COMPLEX*16 array, dimension (1+(N-1)*INCC) The sines of the plane rotations.INCCINCC is INTEGER The increment between elements of C and S. INCC > 0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlarcm(integerM,integerN,doubleprecision,dimension(lda,*)A,integerLDA,complex*16,dimension(ldb,*)B,integerLDB,complex*16,dimension(ldc,*)C,integerLDC,doubleprecision,dimension(*)RWORK)ZLARCMcopies all or part of a real two-dimensional array to a complex array.Purpose:ZLARCM performs a very simple matrix-matrix multiplication: C := A * B, where A is M by M and real; B is M by N and complex; C is M by N and complex.Parameters:MM is INTEGER The number of rows of the matrix A and of the matrix C. M >= 0.NN is INTEGER The number of columns and rows of the matrix B and the number of columns of the matrix C. N >= 0.AA is DOUBLE PRECISION array, dimension (LDA, M) On entry, A contains the M by M matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >=max(1,M).BB is COMPLEX*16 array, dimension (LDB, N) On entry, B contains the M by N matrix B.LDBLDB is INTEGER The leading dimension of the array B. LDB >=max(1,M).CC is COMPLEX*16 array, dimension (LDC, N) On exit, C contains the M by N matrix C.LDCLDC is INTEGER The leading dimension of the array C. LDC >=max(1,M).RWORKRWORK is DOUBLE PRECISION array, dimension (2*M*N)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinezlarf(characterSIDE,integerM,integerN,complex*16,dimension(*)V,integerINCV,complex*16TAU,complex*16,dimension(ldc,*)C,integerLDC,complex*16,dimension(*)WORK)ZLARFapplies an elementary reflector to a general rectangular matrix.Purpose:ZLARF applies a complex elementary reflector H to a complex M-by-N matrix C, from either the left or the right. H is represented in the form H = I - tau * v * v**H where tau is a complex scalar and v is a complex vector. If tau = 0, then H is taken to be the unit matrix. To apply H**H, supply conjg(tau) instead tau.Parameters:SIDESIDE is CHARACTER*1 = 'L': form H * C = 'R': form C * HMM is INTEGER The number of rows of the matrix C.NN is INTEGER The number of columns of the matrix C.VV is COMPLEX*16 array, dimension (1 + (M-1)*abs(INCV)) if SIDE = 'L' or (1 + (N-1)*abs(INCV)) if SIDE = 'R' The vector v in the representation of H. V is not used if TAU = 0.INCVINCV is INTEGER The increment between elements of v. INCV <> 0.TAUTAU is COMPLEX*16 The value tau in the representation of H.CC 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'.LDCLDC is INTEGER The leading dimension of the array C. LDC >= max(1,M).WORKWORK is COMPLEX*16 array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R'Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlarfb(characterSIDE,characterTRANS,characterDIRECT,characterSTOREV,integerM,integerN,integerK,complex*16,dimension(ldv,*)V,integerLDV,complex*16,dimension(ldt,*)T,integerLDT,complex*16,dimension(ldc,*)C,integerLDC,complex*16,dimension(ldwork,*)WORK,integerLDWORK)ZLARFBapplies a block reflector or its conjugate-transpose to a general rectangular matrix.Purpose:ZLARFB applies a complex block reflector H or its transpose H**H to a complex M-by-N matrix C, from either the left or the right.Parameters:SIDESIDE is CHARACTER*1 = 'L': apply H or H**H from the Left = 'R': apply H or H**H from the RightTRANSTRANS is CHARACTER*1 = 'N': apply H (No transpose) = 'C': apply H**H (Conjugate transpose)DIRECTDIRECT is CHARACTER*1 Indicates how H is formed from a product of elementary reflectors = 'F': H = H(1) H(2) . . . H(k) (Forward) = 'B': H = H(k) . . . H(2) H(1) (Backward)STOREVSTOREV is CHARACTER*1 Indicates how the vectors which define the elementary reflectors are stored: = 'C': Columnwise = 'R': RowwiseMM is INTEGER The number of rows of the matrix C.NN is INTEGER The number of columns of the matrix C.KK is INTEGER The order of the matrix T (= the number of elementary reflectors whose product defines the block reflector).VV is COMPLEX*16 array, dimension (LDV,K) if STOREV = 'C' (LDV,M) if STOREV = 'R' and SIDE = 'L' (LDV,N) if STOREV = 'R' and SIDE = 'R' See Further Details.LDVLDV is INTEGER The leading dimension of the array V. If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); if STOREV = 'R', LDV >= K.TT is COMPLEX*16 array, dimension (LDT,K) The triangular K-by-K matrix T in the representation of the block reflector.LDTLDT is INTEGER The leading dimension of the array T. LDT >= K.CC 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.LDCLDC is INTEGER The leading dimension of the array C. LDC >= max(1,M).WORKWORK is COMPLEX*16 array, dimension (LDWORK,K)LDWORKLDWORK is INTEGER The leading dimension of the array WORK. If SIDE = 'L', LDWORK >= max(1,N); if SIDE = 'R', LDWORK >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2013FurtherDetails:The shape of the matrix V and the storage of the vectors which define the H(i) is best illustrated by the following example with n = 5 and k = 3. The elements equal to 1 are not stored; the corresponding array elements are modified but restored on exit. The rest of the array is not used. DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) ( v1 1 ) ( 1 v2 v2 v2 ) ( v1 v2 1 ) ( 1 v3 v3 ) ( v1 v2 v3 ) ( v1 v2 v3 ) DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': V = ( v1 v2 v3 ) V = ( v1 v1 1 ) ( v1 v2 v3 ) ( v2 v2 v2 1 ) ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) ( 1 v3 ) ( 1 )subroutinezlarfg(integerN,complex*16ALPHA,complex*16,dimension(*)X,integerINCX,complex*16TAU)ZLARFGgenerates an elementary reflector (Householder matrix).Purpose:ZLARFG generates a complex elementary reflector H of order n, such that H**H * ( alpha ) = ( beta ), H**H * H = I. ( x ) ( 0 ) where alpha and beta are scalars, with beta real, and x is an (n-1)-element complex vector. H is represented in the form H = I - tau * ( 1 ) * ( 1 v**H ) , ( v ) where tau is a complex scalar and v is a complex (n-1)-element vector. Note that H is not hermitian. If the elements of x are all zero and alpha is real, then tau = 0 and H is taken to be the unit matrix. Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 .Parameters:NN is INTEGER The order of the elementary reflector.ALPHAALPHA is COMPLEX*16 On entry, the value alpha. On exit, it is overwritten with the value beta.XX is COMPLEX*16 array, dimension (1+(N-2)*abs(INCX)) On entry, the vector x. On exit, it is overwritten with the vector v.INCXINCX is INTEGER The increment between elements of X. INCX > 0.TAUTAU is COMPLEX*16 The value tau.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:November 2017subroutinezlarfgp(integerN,complex*16ALPHA,complex*16,dimension(*)X,integerINCX,complex*16TAU)ZLARFGPgenerates an elementary reflector (Householder matrix) with non-negative beta.Purpose:ZLARFGP generates a complex elementary reflector H of order n, such that H**H * ( alpha ) = ( beta ), H**H * H = I. ( x ) ( 0 ) where alpha and beta are scalars, beta is real and non-negative, and x is an (n-1)-element complex vector. H is represented in the form H = I - tau * ( 1 ) * ( 1 v**H ) , ( v ) where tau is a complex scalar and v is a complex (n-1)-element vector. Note that H is not hermitian. If the elements of x are all zero and alpha is real, then tau = 0 and H is taken to be the unit matrix.Parameters:NN is INTEGER The order of the elementary reflector.ALPHAALPHA is COMPLEX*16 On entry, the value alpha. On exit, it is overwritten with the value beta.XX is COMPLEX*16 array, dimension (1+(N-2)*abs(INCX)) On entry, the vector x. On exit, it is overwritten with the vector v.INCXINCX is INTEGER The increment between elements of X. INCX > 0.TAUTAU is COMPLEX*16 The value tau.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:November 2017subroutinezlarft(characterDIRECT,characterSTOREV,integerN,integerK,complex*16,dimension(ldv,*)V,integerLDV,complex*16,dimension(*)TAU,complex*16,dimension(ldt,*)T,integerLDT)ZLARFTforms the triangular factor T of a block reflector H = I - vtvHPurpose:ZLARFT forms the triangular factor T of a complex block reflector H of order n, which is defined as a product of k elementary reflectors. If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. If STOREV = 'C', the vector which defines the elementary reflector H(i) is stored in the i-th column of the array V, and H = I - V * T * V**H If STOREV = 'R', the vector which defines the elementary reflector H(i) is stored in the i-th row of the array V, and H = I - V**H * T * VParameters:DIRECTDIRECT is CHARACTER*1 Specifies the order in which the elementary reflectors are multiplied to form the block reflector: = 'F': H = H(1) H(2) . . . H(k) (Forward) = 'B': H = H(k) . . . H(2) H(1) (Backward)STOREVSTOREV is CHARACTER*1 Specifies how the vectors which define the elementary reflectors are stored (see also Further Details): = 'C': columnwise = 'R': rowwiseNN is INTEGER The order of the block reflector H. N >= 0.KK is INTEGER The order of the triangular factor T (= the number of elementary reflectors). K >= 1.VV is COMPLEX*16 array, dimension (LDV,K) if STOREV = 'C' (LDV,N) if STOREV = 'R' The matrix V. See further details.LDVLDV is INTEGER The leading dimension of the array V. If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.TAUTAU is COMPLEX*16 array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i).TT 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.LDTLDT is INTEGER The leading dimension of the array T. LDT >= K.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016FurtherDetails:The shape of the matrix V and the storage of the vectors which define the H(i) is best illustrated by the following example with n = 5 and k = 3. The elements equal to 1 are not stored. DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) ( v1 1 ) ( 1 v2 v2 v2 ) ( v1 v2 1 ) ( 1 v3 v3 ) ( v1 v2 v3 ) ( v1 v2 v3 ) DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': V = ( v1 v2 v3 ) V = ( v1 v1 1 ) ( v1 v2 v3 ) ( v2 v2 v2 1 ) ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) ( 1 v3 ) ( 1 )subroutinezlarfx(characterSIDE,integerM,integerN,complex*16,dimension(*)V,complex*16TAU,complex*16,dimension(ldc,*)C,integerLDC,complex*16,dimension(*)WORK)ZLARFXapplies an elementary reflector to a general rectangular matrix, with loop unrolling when the reflector has order ≤ 10.Purpose:ZLARFX applies a complex elementary reflector H to a complex m by n matrix C, from either the left or the right. H is represented in the form H = I - tau * v * v**H where tau is a complex scalar and v is a complex vector. If tau = 0, then H is taken to be the unit matrix This version uses inline code if H has order < 11.Parameters:SIDESIDE is CHARACTER*1 = 'L': form H * C = 'R': form C * HMM is INTEGER The number of rows of the matrix C.NN is INTEGER The number of columns of the matrix C.VV is COMPLEX*16 array, dimension (M) if SIDE = 'L' or (N) if SIDE = 'R' The vector v in the representation of H.TAUTAU is COMPLEX*16 The value tau in the representation of H.CC 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'.LDCLDC is INTEGER The leading dimension of the array C. LDA >= max(1,M).WORKWORK is COMPLEX*16 array, dimension (N) if SIDE = 'L' or (M) if SIDE = 'R' WORK is not referenced if H has order < 11.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlargv(integerN,complex*16,dimension(*)X,integerINCX,complex*16,dimension(*)Y,integerINCY,doubleprecision,dimension(*)C,integerINCC)ZLARGVgenerates a vector of plane rotations with real cosines and complex sines.Purpose:ZLARGV generates a vector of complex plane rotations with real cosines, determined by elements of the complex vectors x and y. For i = 1,2,...,n ( c(i) s(i) ) ( x(i) ) = ( r(i) ) ( -conjg(s(i)) c(i) ) ( y(i) ) = ( 0 ) where c(i)**2 + ABS(s(i))**2 = 1 The following conventions are used (these are the same as in ZLARTG, but differ from the BLAS1 routine ZROTG): If y(i)=0, then c(i)=1 and s(i)=0. If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.Parameters:NN is INTEGER The number of plane rotations to be generated.XX is COMPLEX*16 array, dimension (1+(N-1)*INCX) On entry, the vector x. On exit, x(i) is overwritten by r(i), for i = 1,...,n.INCXINCX is INTEGER The increment between elements of X. INCX > 0.YY is COMPLEX*16 array, dimension (1+(N-1)*INCY) On entry, the vector y. On exit, the sines of the plane rotations.INCYINCY is INTEGER The increment between elements of Y. INCY > 0.CC is DOUBLE PRECISION array, dimension (1+(N-1)*INCC) The cosines of the plane rotations.INCCINCC is INTEGER The increment between elements of C. INCC > 0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel This version has a few statements commented out for thread safety (machine parameters are computed on each entry). 10 feb 03, SJH.subroutinezlarnv(integerIDIST,integer,dimension(4)ISEED,integerN,complex*16,dimension(*)X)ZLARNVreturns a vector of random numbers from a uniform or normal distribution.Purpose:ZLARNV returns a vector of n random complex numbers from a uniform or normal distribution.Parameters:IDISTIDIST is INTEGER Specifies the distribution of the random numbers: = 1: real and imaginary parts each uniform (0,1) = 2: real and imaginary parts each uniform (-1,1) = 3: real and imaginary parts each normal (0,1) = 4: uniformly distributed on the disc abs(z) < 1 = 5: uniformly distributed on the circle abs(z) = 1ISEEDISEED is INTEGER array, dimension (4) On entry, the seed of the random number generator; the array elements must be between 0 and 4095, and ISEED(4) must be odd. On exit, the seed is updated.NN is INTEGER The number of random numbers to be generated.XX is COMPLEX*16 array, dimension (N) The generated random numbers.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:This routine calls the auxiliary routine DLARUV to generate random real numbers from a uniform (0,1) distribution, in batches of up to 128 using vectorisable code. The Box-Muller method is used to transform numbers from a uniform to a normal distribution.subroutinezlarrv(integerN,doubleprecisionVL,doubleprecisionVU,doubleprecision,dimension(*)D,doubleprecision,dimension(*)L,doubleprecisionPIVMIN,integer,dimension(*)ISPLIT,integerM,integerDOL,integerDOU,doubleprecisionMINRGP,doubleprecisionRTOL1,doubleprecisionRTOL2,doubleprecision,dimension(*)W,doubleprecision,dimension(*)WERR,doubleprecision,dimension(*)WGAP,integer,dimension(*)IBLOCK,integer,dimension(*)INDEXW,doubleprecision,dimension(*)GERS,complex*16,dimension(ldz,*)Z,integerLDZ,integer,dimension(*)ISUPPZ,doubleprecision,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)ZLARRVcomputes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT.Purpose:ZLARRV computes the eigenvectors of the tridiagonal matrix T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. The input eigenvalues should have been computed by DLARRE.Parameters:NN is INTEGER The order of the matrix. N >= 0.VLVL is DOUBLE PRECISION Lower bound of the interval that contains the desired eigenvalues. VL < VU. Needed to compute gaps on the left or right end of the extremal eigenvalues in the desired RANGE.VUVU is DOUBLE PRECISION Upper bound of the interval that contains the desired eigenvalues. VL < VU. Needed to compute gaps on the left or right end of the extremal eigenvalues in the desired RANGE.DD is DOUBLE PRECISION array, dimension (N) On entry, the N diagonal elements of the diagonal matrix D. On exit, D may be overwritten.LL is DOUBLE PRECISION array, dimension (N) On entry, the (N-1) subdiagonal elements of the unit bidiagonal matrix L are in elements 1 to N-1 of L (if the matrix is not split.) At the end of each block is stored the corresponding shift as given by DLARRE. On exit, L is overwritten.PIVMINPIVMIN is DOUBLE PRECISION The minimum pivot allowed in the Sturm sequence.ISPLITISPLIT is INTEGER array, dimension (N) The splitting points, at which T breaks up into blocks. The first block consists of rows/columns 1 to ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 through ISPLIT( 2 ), etc.MM is INTEGER The total number of input eigenvalues. 0 <= M <= N.DOLDOL is INTEGERDOUDOU is INTEGER If the user wants to compute only selected eigenvectors from all the eigenvalues supplied, he can specify an index range DOL:DOU. Or else the setting DOL=1, DOU=M should be applied. Note that DOL and DOU refer to the order in which the eigenvalues are stored in W. If the user wants to compute only selected eigenpairs, then the columns DOL-1 to DOU+1 of the eigenvector space Z contain the computed eigenvectors. All other columns of Z are set to zero.MINRGPMINRGP is DOUBLE PRECISIONRTOL1RTOL1 is DOUBLE PRECISIONRTOL2RTOL2 is DOUBLE PRECISION Parameters for bisection. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )WW is DOUBLE PRECISION array, dimension (N) The first M elements of W contain the APPROXIMATE eigenvalues for which eigenvectors are to be computed. The eigenvalues should be grouped by split-off block and ordered from smallest to largest within the block ( The output array W from DLARRE is expected here ). Furthermore, they are with respect to the shift of the corresponding root representation for their block. On exit, W holds the eigenvalues of the UNshifted matrix.WERRWERR is DOUBLE PRECISION array, dimension (N) The first M elements contain the semiwidth of the uncertainty interval of the corresponding eigenvalue in WWGAPWGAP is DOUBLE PRECISION array, dimension (N) The separation from the right neighbor eigenvalue in W.IBLOCKIBLOCK is INTEGER array, dimension (N) The indices of the blocks (submatrices) associated with the corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to the first block from the top, =2 if W(i) belongs to the second block, etc.INDEXWINDEXW is INTEGER array, dimension (N) The indices of the eigenvalues within each block (submatrix); for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.GERSGERS is DOUBLE PRECISION array, dimension (2*N) The N Gerschgorin intervals (the i-th Gerschgorin interval is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should be computed from the original UNshifted matrix.ZZ is COMPLEX*16 array, dimension (LDZ, max(1,M) ) If INFO = 0, the first M columns of Z contain the orthonormal eigenvectors of the matrix T corresponding to the input eigenvalues, with the i-th column of Z holding the eigenvector associated with W(i). Note: the user must ensure that at least max(1,M) columns are supplied in the array Z.LDZLDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', LDZ >= max(1,N).ISUPPZISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) The support of the eigenvectors in Z, i.e., the indices indicating the nonzero elements in Z. The I-th eigenvector is nonzero only in elements ISUPPZ( 2*I-1 ) through ISUPPZ( 2*I ).WORKWORK is DOUBLE PRECISION array, dimension (12*N)IWORKIWORK is INTEGER array, dimension (7*N)INFOINFO is INTEGER = 0: successful exit > 0: A problem occurred in ZLARRV. < 0: One of the called subroutines signaled an internal problem. Needs inspection of the corresponding parameter IINFO for further information. =-1: Problem in DLARRB when refining a child's eigenvalues. =-2: Problem in DLARRF when computing the RRR of a child. When a child is inside a tight cluster, it can be difficult to find an RRR. A partial remedy from the user's point of view is to make the parameter MINRGP smaller and recompile. However, as the orthogonality of the computed vectors is proportional to 1/MINRGP, the user should be aware that he might be trading in precision when he decreases MINRGP. =-3: Problem in DLARRB when refining a single eigenvalue after the Rayleigh correction was rejected. = 5: The Rayleigh Quotient Iteration failed to converge to full accuracy in MAXITR steps.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Beresford Parlett, University of California, Berkeley, USA Jim Demmel, University of California, Berkeley, USA Inderjit Dhillon, University of Texas, Austin, USA Osni Marques, LBNL/NERSC, USA Christof Voemel, University of California, Berkeley, USAsubroutinezlartg(complex*16F,complex*16G,doubleprecisionCS,complex*16SN,complex*16R)ZLARTGgenerates a plane rotation with real cosine and complex sine.Purpose:ZLARTG generates a plane rotation so that [ CS SN ] [ F ] [ R ] [ __ ] . [ ] = [ ] where CS**2 + |SN|**2 = 1. [ -SN CS ] [ G ] [ 0 ] This is a faster version of the BLAS1 routine ZROTG, except for the following differences: F and G are unchanged on return. If G=0, then CS=1 and SN=0. If F=0, then CS=0 and SN is chosen so that R is real.Parameters:FF is COMPLEX*16 The first component of vector to be rotated.GG is COMPLEX*16 The second component of vector to be rotated.CSCS is DOUBLE PRECISION The cosine of the rotation.SNSN is COMPLEX*16 The sine of the rotation.RR is COMPLEX*16 The nonzero component of the rotated vector.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel This version has a few statements commented out for thread safety (machine parameters are computed on each entry). 10 feb 03, SJH.subroutinezlartv(integerN,complex*16,dimension(*)X,integerINCX,complex*16,dimension(*)Y,integerINCY,doubleprecision,dimension(*)C,complex*16,dimension(*)S,integerINCC)ZLARTVapplies a vector of plane rotations with real cosines and complex sines to the elements of a pair of vectors.Purpose:ZLARTV applies a vector of complex plane rotations with real cosines to elements of the complex vectors x and y. For i = 1,2,...,n ( x(i) ) := ( c(i) s(i) ) ( x(i) ) ( y(i) ) ( -conjg(s(i)) c(i) ) ( y(i) )Parameters:NN is INTEGER The number of plane rotations to be applied.XX is COMPLEX*16 array, dimension (1+(N-1)*INCX) The vector x.INCXINCX is INTEGER The increment between elements of X. INCX > 0.YY is COMPLEX*16 array, dimension (1+(N-1)*INCY) The vector y.INCYINCY is INTEGER The increment between elements of Y. INCY > 0.CC is DOUBLE PRECISION array, dimension (1+(N-1)*INCC) The cosines of the plane rotations.SS is COMPLEX*16 array, dimension (1+(N-1)*INCC) The sines of the plane rotations.INCCINCC is INTEGER The increment between elements of C and S. INCC > 0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlascl(characterTYPE,integerKL,integerKU,doubleprecisionCFROM,doubleprecisionCTO,integerM,integerN,complex*16,dimension(lda,*)A,integerLDA,integerINFO)ZLASCLmultiplies a general rectangular matrix by a real scalar defined as cto/cfrom.Purpose:ZLASCL multiplies the M by N complex matrix A by the real scalar CTO/CFROM. This is done without over/underflow as long as the final result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that A may be full, upper triangular, lower triangular, upper Hessenberg, or banded.Parameters:TYPETYPE is CHARACTER*1 TYPE indices the storage type of the input matrix. = 'G': A is a full matrix. = 'L': A is a lower triangular matrix. = 'U': A is an upper triangular matrix. = 'H': A is an upper Hessenberg matrix. = 'B': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the lower half stored. = 'Q': A is a symmetric band matrix with lower bandwidth KL and upper bandwidth KU and with the only the upper half stored. = 'Z': A is a band matrix with lower bandwidth KL and upper bandwidth KU. See ZGBTRF for storage details.KLKL is INTEGER The lower bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'.KUKU is INTEGER The upper bandwidth of A. Referenced only if TYPE = 'B', 'Q' or 'Z'.CFROMCFROM is DOUBLE PRECISIONCTOCTO is DOUBLE PRECISION The matrix A is multiplied by CTO/CFROM. A(I,J) is computed without over/underflow if the final result CTO*A(I,J)/CFROM can be represented without over/underflow. CFROM must be nonzero.MM is INTEGER The number of rows of the matrix A. M >= 0.NN is INTEGER The number of columns of the matrix A. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) The matrix to be multiplied by CTO/CFROM. See TYPE for the storage type.LDALDA is INTEGER The leading dimension of the array A. If TYPE = 'G', 'L', 'U', 'H', LDA >= max(1,M); TYPE = 'B', LDA >= KL+1; TYPE = 'Q', LDA >= KU+1; TYPE = 'Z', LDA >= 2*KL+KU+1.INFOINFO is INTEGER 0 - successful exit <0 - if INFO = -i, the i-th argument had an illegal value.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinezlaset(characterUPLO,integerM,integerN,complex*16ALPHA,complex*16BETA,complex*16,dimension(lda,*)A,integerLDA)ZLASETinitializes the off-diagonal elements and the diagonal elements of a matrix to given values.Purpose:ZLASET initializes a 2-D array A to BETA on the diagonal and ALPHA on the offdiagonals.Parameters:UPLOUPLO is CHARACTER*1 Specifies the part of the matrix A to be set. = 'U': Upper triangular part is set. The lower triangle is unchanged. = 'L': Lower triangular part is set. The upper triangle is unchanged. Otherwise: All of the matrix A is set.MM is INTEGER On entry, M specifies the number of rows of A.NN is INTEGER On entry, N specifies the number of columns of A.ALPHAALPHA is COMPLEX*16 All the offdiagonal array elements are set to ALPHA.BETABETA is COMPLEX*16 All the diagonal array elements are set to BETA.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the m by n matrix A. On exit, A(i,j) = ALPHA, 1 <= i <= m, 1 <= j <= n, i.ne.j; A(i,i) = BETA , 1 <= i <= min(m,n)LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlasr(characterSIDE,characterPIVOT,characterDIRECT,integerM,integerN,doubleprecision,dimension(*)C,doubleprecision,dimension(*)S,complex*16,dimension(lda,*)A,integerLDA)ZLASRapplies a sequence of plane rotations to a general rectangular matrix.Purpose:ZLASR applies a sequence of real plane rotations to a complex matrix A, from either the left or the right. When SIDE = 'L', the transformation takes the form A := P*A and when SIDE = 'R', the transformation takes the form A := A*P**T where P is an orthogonal matrix consisting of a sequence of z plane rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', and P**T is the transpose of P. When DIRECT = 'F' (Forward sequence), then P = P(z-1) * ... * P(2) * P(1) and when DIRECT = 'B' (Backward sequence), then P = P(1) * P(2) * ... * P(z-1) where P(k) is a plane rotation matrix defined by the 2-by-2 rotation R(k) = ( c(k) s(k) ) = ( -s(k) c(k) ). When PIVOT = 'V' (Variable pivot), the rotation is performed for the plane (k,k+1), i.e., P(k) has the form P(k) = ( 1 ) ( ... ) ( 1 ) ( c(k) s(k) ) ( -s(k) c(k) ) ( 1 ) ( ... ) ( 1 ) where R(k) appears as a rank-2 modification to the identity matrix in rows and columns k and k+1. When PIVOT = 'T' (Top pivot), the rotation is performed for the plane (1,k+1), so P(k) has the form P(k) = ( c(k) s(k) ) ( 1 ) ( ... ) ( 1 ) ( -s(k) c(k) ) ( 1 ) ( ... ) ( 1 ) where R(k) appears in rows and columns 1 and k+1. Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is performed for the plane (k,z), giving P(k) the form P(k) = ( 1 ) ( ... ) ( 1 ) ( c(k) s(k) ) ( 1 ) ( ... ) ( 1 ) ( -s(k) c(k) ) where R(k) appears in rows and columns k and z. The rotations are performed without ever forming P(k) explicitly.Parameters:SIDESIDE is CHARACTER*1 Specifies whether the plane rotation matrix P is applied to A on the left or the right. = 'L': Left, compute A := P*A = 'R': Right, compute A:= A*P**TPIVOTPIVOT is CHARACTER*1 Specifies the plane for which P(k) is a plane rotation matrix. = 'V': Variable pivot, the plane (k,k+1) = 'T': Top pivot, the plane (1,k+1) = 'B': Bottom pivot, the plane (k,z)DIRECTDIRECT is CHARACTER*1 Specifies whether P is a forward or backward sequence of plane rotations. = 'F': Forward, P = P(z-1)*...*P(2)*P(1) = 'B': Backward, P = P(1)*P(2)*...*P(z-1)MM is INTEGER The number of rows of the matrix A. If m <= 1, an immediate return is effected.NN is INTEGER The number of columns of the matrix A. If n <= 1, an immediate return is effected.CC is DOUBLE PRECISION array, dimension (M-1) if SIDE = 'L' (N-1) if SIDE = 'R' The cosines c(k) of the plane rotations.SS is DOUBLE PRECISION array, dimension (M-1) if SIDE = 'L' (N-1) if SIDE = 'R' The sines s(k) of the plane rotations. The 2-by-2 plane rotation part of the matrix P(k), R(k), has the form R(k) = ( c(k) s(k) ) ( -s(k) c(k) ).AA is COMPLEX*16 array, dimension (LDA,N) The M-by-N matrix A. On exit, A is overwritten by P*A if SIDE = 'R' or by A*P**T if SIDE = 'L'.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlassq(integerN,complex*16,dimension(*)X,integerINCX,doubleprecisionSCALE,doubleprecisionSUMSQ)ZLASSQupdates a sum of squares represented in scaled form.Purpose:ZLASSQ returns the values scl and ssq such that ( scl**2 )*ssq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is assumed to be at least unity and the value of ssq will then satisfy 1.0 .le. ssq .le. ( sumsq + 2*n ). scale is assumed to be non-negative and scl returns the value scl = max( scale, abs( real( x( i ) ) ), abs( aimag( x( i ) ) ) ), i scale and sumsq must be supplied in SCALE and SUMSQ respectively. SCALE and SUMSQ are overwritten by scl and ssq respectively. The routine makes only one pass through the vector X.Parameters:NN is INTEGER The number of elements to be used from the vector X.XX is COMPLEX*16 array, dimension (N) The vector x as described above. x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.INCXINCX is INTEGER The increment between successive values of the vector X. INCX > 0.SCALESCALE is DOUBLE PRECISION On entry, the value scale in the equation above. On exit, SCALE is overwritten with the value scl .SUMSQSUMSQ is DOUBLE PRECISION On entry, the value sumsq in the equation above. On exit, SUMSQ is overwritten with the value ssq .Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlaswp(integerN,complex*16,dimension(lda,*)A,integerLDA,integerK1,integerK2,integer,dimension(*)IPIV,integerINCX)ZLASWPperforms a series of row interchanges on a general rectangular matrix.Purpose:ZLASWP performs a series of row interchanges on the matrix A. One row interchange is initiated for each of rows K1 through K2 of A.Parameters:NN is INTEGER The number of columns of the matrix A.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the matrix of column dimension N to which the row interchanges will be applied. On exit, the permuted matrix.LDALDA is INTEGER The leading dimension of the array A.K1K1 is INTEGER The first element of IPIV for which a row interchange will be done.K2K2 is INTEGER (K2-K1+1) is the number of elements of IPIV for which a row interchange will be done.IPIVIPIV is INTEGER array, dimension (K1+(K2-K1)*abs(INCX)) The vector of pivot indices. Only the elements in positions K1 through K1+(K2-K1)*abs(INCX) of IPIV are accessed. IPIV(K1+(K-K1)*abs(INCX)) = L implies rows K and L are to be interchanged.INCXINCX is INTEGER The increment between successive values of IPIV. If INCX is negative, the pivots are applied in reverse order.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017FurtherDetails:Modified by R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USAsubroutinezlat2c(characterUPLO,integerN,complex*16,dimension(lda,*)A,integerLDA,complex,dimension(ldsa,*)SA,integerLDSA,integerINFO)ZLAT2Cconverts a double complex triangular matrix to a complex triangular matrix.Purpose:ZLAT2C converts a COMPLEX*16 triangular matrix, SA, to a COMPLEX triangular matrix, A. RMAX is the overflow for the SINGLE PRECISION arithmetic ZLAT2C checks that all the entries of A are between -RMAX and RMAX. If not the conversion is aborted and a flag is raised. This is an auxiliary routine so there is no argument checking.Parameters:UPLOUPLO is CHARACTER*1 = 'U': A is upper triangular; = 'L': A is lower triangular.NN is INTEGER The number of rows and columns of the matrix A. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the N-by-N triangular coefficient matrix A.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).SASA is COMPLEX array, dimension (LDSA,N) Only the UPLO part of SA is referenced. On exit, if INFO=0, the N-by-N coefficient matrix SA; if INFO>0, the content of the UPLO part of SA is unspecified.LDSALDSA is INTEGER The leading dimension of the array SA. LDSA >= max(1,M).INFOINFO is INTEGER = 0: successful exit. = 1: an entry of the matrix A is greater than the SINGLE PRECISION overflow threshold, in this case, the content of the UPLO part of SA in exit is unspecified.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlatbs(characterUPLO,characterTRANS,characterDIAG,characterNORMIN,integerN,integerKD,complex*16,dimension(ldab,*)AB,integerLDAB,complex*16,dimension(*)X,doubleprecisionSCALE,doubleprecision,dimension(*)CNORM,integerINFO)ZLATBSsolves a triangular banded system of equations.Purpose:ZLATBS solves one of the triangular systems A * x = s*b, A**T * x = s*b, or A**H * x = s*b, with scaling to prevent overflow, where A is an upper or lower triangular band matrix. Here A**T denotes the transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine ZTBSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangularTRANSTRANS is CHARACTER*1 Specifies the operation applied to A. = 'N': Solve A * x = s*b (No transpose) = 'T': Solve A**T * x = s*b (Transpose) = 'C': Solve A**H * x = s*b (Conjugate transpose)DIAGDIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangularNORMINNORMIN is CHARACTER*1 Specifies whether CNORM has been set or not. = 'Y': CNORM contains the column norms on entry = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.NN is INTEGER The order of the matrix A. N >= 0.KDKD is INTEGER The number of subdiagonals or superdiagonals in the triangular matrix A. KD >= 0.ABAB is COMPLEX*16 array, dimension (LDAB,N) The upper or lower triangular band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).LDABLDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.XX is COMPLEX*16 array, dimension (N) On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.SCALESCALE is DOUBLE PRECISION The scaling factor s for the triangular system A * x = s*b, A**T * x = s*b, or A**H * x = s*b. If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.CNORMCNORM is DOUBLE PRECISION array, dimension (N) If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm. If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:November 2017FurtherDetails:A rough bound on x is computed; if that is less than overflow, ZTBSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation. A columnwise scheme is used for solving A*x = b. The basic algorithm if A is lower triangular is x[1:n] := b[1:n] for j = 1, ..., n x(j) := x(j) / A(j,j) x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] end Define bounds on the components of x after j iterations of the loop: M(j) = bound on x[1:j] G(j) = bound on x[j+1:n] Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. Then for iteration j+1 we have M(j+1) <= G(j) / | A(j+1,j+1) | G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) where CNORM(j+1) is greater than or equal to the infinity-norm of column j+1 of A, not counting the diagonal. Hence G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) 1<=i<=j and |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) 1<=i< j Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTBSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow). The bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. If the computed bound is greater than a large constant, x is scaled to prevent overflow, but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *x = b. The basic algorithm for A upper triangular is for j = 1, ..., n x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) end We simultaneously compute two bounds G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j M(j) = bound on x(i), 1<=i<=j The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. Then the bound on x(j) is M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 1<=i<=j and we can safely call ZTBSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).subroutinezlatdf(integerIJOB,integerN,complex*16,dimension(ldz,*)Z,integerLDZ,complex*16,dimension(*)RHS,doubleprecisionRDSUM,doubleprecisionRDSCAL,integer,dimension(*)IPIV,integer,dimension(*)JPIV)ZLATDFuses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate.Purpose:ZLATDF computes the contribution to the reciprocal Dif-estimate by solving for x in Z * x = b, where b is chosen such that the norm of x is as large as possible. It is assumed that LU decomposition of Z has been computed by ZGETC2. On entry RHS = f holds the contribution from earlier solved sub-systems, and on return RHS = x. The factorization of Z returned by ZGETC2 has the form Z = P * L * U * Q, where P and Q are permutation matrices. L is lower triangular with unit diagonal elements and U is upper triangular.Parameters:IJOBIJOB is INTEGER IJOB = 2: First compute an approximative null-vector e of Z using ZGECON, e is normalized and solve for Zx = +-e - f with the sign giving the greater value of 2-norm(x). About 5 times as expensive as Default. IJOB .ne. 2: Local look ahead strategy where all entries of the r.h.s. b is chosen as either +1 or -1. Default.NN is INTEGER The number of columns of the matrix Z.ZZ is COMPLEX*16 array, dimension (LDZ, N) On entry, the LU part of the factorization of the n-by-n matrix Z computed by ZGETC2: Z = P * L * U * QLDZLDZ is INTEGER The leading dimension of the array Z. LDA >= max(1, N).RHSRHS is COMPLEX*16 array, dimension (N). On entry, RHS contains contributions from other subsystems. On exit, RHS contains the solution of the subsystem with entries according to the value of IJOB (see above).RDSUMRDSUM is DOUBLE PRECISION On entry, the sum of squares of computed contributions to the Dif-estimate under computation by ZTGSYL, where the scaling factor RDSCAL (see below) has been factored out. On exit, the corresponding sum of squares updated with the contributions from the current sub-system. If TRANS = 'T' RDSUM is not touched. NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL.RDSCALRDSCAL is DOUBLE PRECISION On entry, scaling factor used to prevent overflow in RDSUM. On exit, RDSCAL is updated w.r.t. the current contributions in RDSUM. If TRANS = 'T', RDSCAL is not touched. NOTE: RDSCAL only makes sense when ZTGSY2 is called by ZTGSYL.IPIVIPIV is INTEGER array, dimension (N). The pivot indices; for 1 <= i <= N, row i of the matrix has been interchanged with row IPIV(i).JPIVJPIV is INTEGER array, dimension (N). The pivot indices; for 1 <= j <= N, column j of the matrix has been interchanged with column JPIV(j).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016FurtherDetails:This routine is a further developed implementation of algorithm BSOLVE in [1] using complete pivoting in the LU factorization.Contributors:Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.References:[1] Bo Kagstrom and Lars Westin, Generalized Schur Methods with Condition Estimators for Solving the Generalized Sylvester Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. [2] Peter Poromaa, On Efficient and Robust Estimators for the Separation between two Regular Matrix Pairs with Applications in Condition Estimation. Report UMINF-95.05, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.subroutinezlatps(characterUPLO,characterTRANS,characterDIAG,characterNORMIN,integerN,complex*16,dimension(*)AP,complex*16,dimension(*)X,doubleprecisionSCALE,doubleprecision,dimension(*)CNORM,integerINFO)ZLATPSsolves a triangular system of equations with the matrix held in packed storage.Purpose:ZLATPS solves one of the triangular systems A * x = s*b, A**T * x = s*b, or A**H * x = s*b, with scaling to prevent overflow, where A is an upper or lower triangular matrix stored in packed form. Here A**T denotes the transpose of A, A**H denotes the conjugate transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine ZTPSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangularTRANSTRANS is CHARACTER*1 Specifies the operation applied to A. = 'N': Solve A * x = s*b (No transpose) = 'T': Solve A**T * x = s*b (Transpose) = 'C': Solve A**H * x = s*b (Conjugate transpose)DIAGDIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangularNORMINNORMIN is CHARACTER*1 Specifies whether CNORM has been set or not. = 'Y': CNORM contains the column norms on entry = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.NN is INTEGER The order of the matrix A. N >= 0.APAP is COMPLEX*16 array, dimension (N*(N+1)/2) The upper or lower triangular matrix A, packed columnwise in a linear array. The j-th column of A is stored in the array AP as follows: if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.XX is COMPLEX*16 array, dimension (N) On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.SCALESCALE is DOUBLE PRECISION The scaling factor s for the triangular system A * x = s*b, A**T * x = s*b, or A**H * x = s*b. If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.CNORMCNORM is DOUBLE PRECISION array, dimension (N) If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm. If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:November 2017FurtherDetails:A rough bound on x is computed; if that is less than overflow, ZTPSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation. A columnwise scheme is used for solving A*x = b. The basic algorithm if A is lower triangular is x[1:n] := b[1:n] for j = 1, ..., n x(j) := x(j) / A(j,j) x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] end Define bounds on the components of x after j iterations of the loop: M(j) = bound on x[1:j] G(j) = bound on x[j+1:n] Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. Then for iteration j+1 we have M(j+1) <= G(j) / | A(j+1,j+1) | G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) where CNORM(j+1) is greater than or equal to the infinity-norm of column j+1 of A, not counting the diagonal. Hence G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) 1<=i<=j and |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) 1<=i< j Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTPSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow). The bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. If the computed bound is greater than a large constant, x is scaled to prevent overflow, but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *x = b. The basic algorithm for A upper triangular is for j = 1, ..., n x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) end We simultaneously compute two bounds G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j M(j) = bound on x(i), 1<=i<=j The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. Then the bound on x(j) is M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 1<=i<=j and we can safely call ZTPSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).subroutinezlatrd(characterUPLO,integerN,integerNB,complex*16,dimension(lda,*)A,integerLDA,doubleprecision,dimension(*)E,complex*16,dimension(*)TAU,complex*16,dimension(ldw,*)W,integerLDW)ZLATRDreduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an unitary similarity transformation.Purpose:ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to Hermitian tridiagonal form by a unitary similarity transformation Q**H * A * Q, and returns the matrices V and W which are needed to apply the transformation to the unreduced part of A. If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a matrix, of which the upper triangle is supplied; if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a matrix, of which the lower triangle is supplied. This is an auxiliary routine called by ZHETRD.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the Hermitian matrix A is stored: = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the matrix A.NBNB is INTEGER The number of rows and columns to be reduced.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the Hermitian matrix A. If UPLO = 'U', the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit: if UPLO = 'U', the last NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements above the diagonal with the array TAU, represent the unitary matrix Q as a product of elementary reflectors; if UPLO = 'L', the first NB columns have been reduced to tridiagonal form, with the diagonal elements overwriting the diagonal elements of A; the elements below the diagonal with the array TAU, represent the unitary matrix Q as a product of elementary reflectors. See Further Details.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).EE is DOUBLE PRECISION array, dimension (N-1) If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal elements of the last NB columns of the reduced matrix; if UPLO = 'L', E(1:nb) contains the subdiagonal elements of the first NB columns of the reduced matrix.TAUTAU is COMPLEX*16 array, dimension (N-1) The scalar factors of the elementary reflectors, stored in TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'. See Further Details.WW is COMPLEX*16 array, dimension (LDW,NB) The n-by-nb matrix W required to update the unreduced part of A.LDWLDW is INTEGER The leading dimension of the array W. LDW >= max(1,N).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:If UPLO = 'U', the matrix Q is represented as a product of elementary reflectors Q = H(n) H(n-1) . . . H(n-nb+1). Each H(i) has the form H(i) = I - tau * v * v**H where tau is a complex scalar, and v is a complex vector with v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i), and tau in TAU(i-1). If UPLO = 'L', the matrix Q is represented as a product of elementary reflectors Q = H(1) H(2) . . . H(nb). Each H(i) has the form H(i) = I - tau * v * v**H where tau is a complex scalar, and v is a complex vector with v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), and tau in TAU(i). The elements of the vectors v together form the n-by-nb matrix V which is needed, with W, to apply the transformation to the unreduced part of the matrix, using a Hermitian rank-2k update of the form: A := A - V*W**H - W*V**H. The contents of A on exit are illustrated by the following examples with n = 5 and nb = 2: if UPLO = 'U': if UPLO = 'L': ( a a a v4 v5 ) ( d ) ( a a v4 v5 ) ( 1 d ) ( a 1 v5 ) ( v1 1 a ) ( d 1 ) ( v1 v2 a a ) ( d ) ( v1 v2 a a a ) where d denotes a diagonal element of the reduced matrix, a denotes an element of the original matrix that is unchanged, and vi denotes an element of the vector defining H(i).subroutinezlatrs(characterUPLO,characterTRANS,characterDIAG,characterNORMIN,integerN,complex*16,dimension(lda,*)A,integerLDA,complex*16,dimension(*)X,doubleprecisionSCALE,doubleprecision,dimension(*)CNORM,integerINFO)ZLATRSsolves a triangular system of equations with the scale factor set to prevent overflow.Purpose:ZLATRS solves one of the triangular systems A * x = s*b, A**T * x = s*b, or A**H * x = s*b, with scaling to prevent overflow. Here A is an upper or lower triangular matrix, A**T denotes the transpose of A, A**H denotes the conjugate transpose of A, x and b are n-element vectors, and s is a scaling factor, usually less than or equal to 1, chosen so that the components of x will be less than the overflow threshold. If the unscaled problem will not cause overflow, the Level 2 BLAS routine ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j), then s is set to 0 and a non-trivial solution to A*x = 0 is returned.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the matrix A is upper or lower triangular. = 'U': Upper triangular = 'L': Lower triangularTRANSTRANS is CHARACTER*1 Specifies the operation applied to A. = 'N': Solve A * x = s*b (No transpose) = 'T': Solve A**T * x = s*b (Transpose) = 'C': Solve A**H * x = s*b (Conjugate transpose)DIAGDIAG is CHARACTER*1 Specifies whether or not the matrix A is unit triangular. = 'N': Non-unit triangular = 'U': Unit triangularNORMINNORMIN is CHARACTER*1 Specifies whether CNORM has been set or not. = 'Y': CNORM contains the column norms on entry = 'N': CNORM is not set on entry. On exit, the norms will be computed and stored in CNORM.NN is INTEGER The order of the matrix A. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) The triangular matrix A. If UPLO = 'U', the leading n by n upper triangular part of the array A contains the upper triangular matrix, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of the array A contains the lower triangular matrix, and the strictly upper triangular part of A is not referenced. If DIAG = 'U', the diagonal elements of A are also not referenced and are assumed to be 1.LDALDA is INTEGER The leading dimension of the array A. LDA >= max (1,N).XX is COMPLEX*16 array, dimension (N) On entry, the right hand side b of the triangular system. On exit, X is overwritten by the solution vector x.SCALESCALE is DOUBLE PRECISION The scaling factor s for the triangular system A * x = s*b, A**T * x = s*b, or A**H * x = s*b. If SCALE = 0, the matrix A is singular or badly scaled, and the vector x is an exact or approximate solution to A*x = 0.CNORMCNORM is DOUBLE PRECISION array, dimension (N) If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains the norm of the off-diagonal part of the j-th column of A. If TRANS = 'N', CNORM(j) must be greater than or equal to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be greater than or equal to the 1-norm. If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns the 1-norm of the offdiagonal part of the j-th column of A.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:November 2017FurtherDetails:A rough bound on x is computed; if that is less than overflow, ZTRSV is called, otherwise, specific code is used which checks for possible overflow or divide-by-zero at every operation. A columnwise scheme is used for solving A*x = b. The basic algorithm if A is lower triangular is x[1:n] := b[1:n] for j = 1, ..., n x(j) := x(j) / A(j,j) x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] end Define bounds on the components of x after j iterations of the loop: M(j) = bound on x[1:j] G(j) = bound on x[j+1:n] Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. Then for iteration j+1 we have M(j+1) <= G(j) / | A(j+1,j+1) | G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) where CNORM(j+1) is greater than or equal to the infinity-norm of column j+1 of A, not counting the diagonal. Hence G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) 1<=i<=j and |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) 1<=i< j Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTRSV if the reciprocal of the largest M(j), j=1,..,n, is larger than max(underflow, 1/overflow). The bound on x(j) is also used to determine when a step in the columnwise method can be performed without fear of overflow. If the computed bound is greater than a large constant, x is scaled to prevent overflow, but if the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. Similarly, a row-wise scheme is used to solve A**T *x = b or A**H *x = b. The basic algorithm for A upper triangular is for j = 1, ..., n x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) end We simultaneously compute two bounds G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j M(j) = bound on x(i), 1<=i<=j The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. Then the bound on x(j) is M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) 1<=i<=j and we can safely call ZTRSV if 1/M(n) and 1/G(n) are both greater than max(underflow, 1/overflow).subroutinezlauu2(characterUPLO,integerN,complex*16,dimension(lda,*)A,integerLDA,integerINFO)ZLAUU2computes the product UUH or LHL, where U and L are upper or lower triangular matrices (unblocked algorithm).Purpose:ZLAUU2 computes the product U * U**H or L**H * L, where the triangular factor U or L is stored in the upper or lower triangular part of the array A. If UPLO = 'U' or 'u' then the upper triangle of the result is stored, overwriting the factor U in A. If UPLO = 'L' or 'l' then the lower triangle of the result is stored, overwriting the factor L in A. This is the unblocked form of the algorithm, calling Level 2 BLAS.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the triangular factor stored in the array A is upper or lower triangular: = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the triangular factor U or L. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the triangular factor U or L. On exit, if UPLO = 'U', the upper triangle of A is overwritten with the upper triangle of the product U * U**H; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**H * L.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezlauum(characterUPLO,integerN,complex*16,dimension(lda,*)A,integerLDA,integerINFO)ZLAUUMcomputes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm).Purpose:ZLAUUM computes the product U * U**H or L**H * L, where the triangular factor U or L is stored in the upper or lower triangular part of the array A. If UPLO = 'U' or 'u' then the upper triangle of the result is stored, overwriting the factor U in A. If UPLO = 'L' or 'l' then the lower triangle of the result is stored, overwriting the factor L in A. This is the blocked form of the algorithm, calling Level 3 BLAS.Parameters:UPLOUPLO is CHARACTER*1 Specifies whether the triangular factor stored in the array A is upper or lower triangular: = 'U': Upper triangular = 'L': Lower triangularNN is INTEGER The order of the triangular factor U or L. N >= 0.AA is COMPLEX*16 array, dimension (LDA,N) On entry, the triangular factor U or L. On exit, if UPLO = 'U', the upper triangle of A is overwritten with the upper triangle of the product U * U**H; if UPLO = 'L', the lower triangle of A is overwritten with the lower triangle of the product L**H * L.LDALDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezrot(integerN,complex*16,dimension(*)CX,integerINCX,complex*16,dimension(*)CY,integerINCY,doubleprecisionC,complex*16S)ZROTapplies a plane rotation with real cosine and complex sine to a pair of complex vectors.Purpose:ZROT applies a plane rotation, where the cos (C) is real and the sin (S) is complex, and the vectors CX and CY are complex.Parameters:NN is INTEGER The number of elements in the vectors CX and CY.CXCX is COMPLEX*16 array, dimension (N) On input, the vector X. On output, CX is overwritten with C*X + S*Y.INCXINCX is INTEGER The increment between successive values of CY. INCX <> 0.CYCY is COMPLEX*16 array, dimension (N) On input, the vector Y. On output, CY is overwritten with -CONJG(S)*X + C*Y.INCYINCY is INTEGER The increment between successive values of CY. INCX <> 0.CC is DOUBLE PRECISIONSS is COMPLEX*16 C and S define a rotation [ C S ] [ -conjg(S) C ] where C*C + S*CONJG(S) = 1.0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezspmv(characterUPLO,integerN,complex*16ALPHA,complex*16,dimension(*)AP,complex*16,dimension(*)X,integerINCX,complex*16BETA,complex*16,dimension(*)Y,integerINCY)ZSPMVcomputes a matrix-vector product for complex vectors using a complex symmetric packed matrixPurpose:ZSPMV performs the matrix-vector operation y := alpha*A*x + beta*y, where alpha and beta are scalars, x and y are n element vectors and A is an n by n symmetric matrix, supplied in packed form.Parameters:UPLOUPLO is CHARACTER*1 On entry, UPLO specifies whether the upper or lower triangular part of the matrix A is supplied in the packed array AP as follows: UPLO = 'U' or 'u' The upper triangular part of A is supplied in AP. UPLO = 'L' or 'l' The lower triangular part of A is supplied in AP. Unchanged on exit.NN is INTEGER On entry, N specifies the order of the matrix A. N must be at least zero. Unchanged on exit.ALPHAALPHA is COMPLEX*16 On entry, ALPHA specifies the scalar alpha. Unchanged on exit.APAP is COMPLEX*16 array, dimension at least ( ( N*( N + 1 ) )/2 ). Before entry, with UPLO = 'U' or 'u', the array AP must contain the upper triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) respectively, and so on. Before entry, with UPLO = 'L' or 'l', the array AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) respectively, and so on. Unchanged on exit.XX is COMPLEX*16 array, dimension at least ( 1 + ( N - 1 )*abs( INCX ) ). Before entry, the incremented array X must contain the N- element vector x. Unchanged on exit.INCXINCX is INTEGER On entry, INCX specifies the increment for the elements of X. INCX must not be zero. Unchanged on exit.BETABETA is COMPLEX*16 On entry, BETA specifies the scalar beta. When BETA is supplied as zero then Y need not be set on input. Unchanged on exit.YY is COMPLEX*16 array, dimension at least ( 1 + ( N - 1 )*abs( INCY ) ). Before entry, the incremented array Y must contain the n element vector y. On exit, Y is overwritten by the updated vector y.INCYINCY is INTEGER On entry, INCY specifies the increment for the elements of Y. INCY must not be zero. Unchanged on exit.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinezspr(characterUPLO,integerN,complex*16ALPHA,complex*16,dimension(*)X,integerINCX,complex*16,dimension(*)AP)ZSPRperforms the symmetrical rank-1 update of a complex symmetric packed matrix.Purpose:ZSPR performs the symmetric rank 1 operation A := alpha*x*x**H + A, where alpha is a complex scalar, x is an n element vector and A is an n by n symmetric matrix, supplied in packed form.Parameters:UPLOUPLO is CHARACTER*1 On entry, UPLO specifies whether the upper or lower triangular part of the matrix A is supplied in the packed array AP as follows: UPLO = 'U' or 'u' The upper triangular part of A is supplied in AP. UPLO = 'L' or 'l' The lower triangular part of A is supplied in AP. Unchanged on exit.NN is INTEGER On entry, N specifies the order of the matrix A. N must be at least zero. Unchanged on exit.ALPHAALPHA is COMPLEX*16 On entry, ALPHA specifies the scalar alpha. Unchanged on exit.XX is COMPLEX*16 array, dimension at least ( 1 + ( N - 1 )*abs( INCX ) ). Before entry, the incremented array X must contain the N- element vector x. Unchanged on exit.INCXINCX is INTEGER On entry, INCX specifies the increment for the elements of X. INCX must not be zero. Unchanged on exit.APAP is COMPLEX*16 array, dimension at least ( ( N*( N + 1 ) )/2 ). Before entry, with UPLO = 'U' or 'u', the array AP must contain the upper triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) respectively, and so on. On exit, the array AP is overwritten by the upper triangular part of the updated matrix. Before entry, with UPLO = 'L' or 'l', the array AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so that AP( 1 ) contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) respectively, and so on. On exit, the array AP is overwritten by the lower triangular part of the updated matrix. Note that the imaginary parts of the diagonal elements need not be set, they are assumed to be zero, and on exit they are set to zero.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutineztprfb(characterSIDE,characterTRANS,characterDIRECT,characterSTOREV,integerM,integerN,integerK,integerL,complex*16,dimension(ldv,*)V,integerLDV,complex*16,dimension(ldt,*)T,integerLDT,complex*16,dimension(lda,*)A,integerLDA,complex*16,dimension(ldb,*)B,integerLDB,complex*16,dimension(ldwork,*)WORK,integerLDWORK)ZTPRFBapplies a real or complex 'triangular-pentagonal' blocked reflector to a real or complex matrix, which is composed of two blocks.Purpose:ZTPRFB applies a complex "triangular-pentagonal" block reflector H or its conjugate transpose H**H to a complex matrix C, which is composed of two blocks A and B, either from the left or right.Parameters:SIDESIDE is CHARACTER*1 = 'L': apply H or H**H from the Left = 'R': apply H or H**H from the RightTRANSTRANS is CHARACTER*1 = 'N': apply H (No transpose) = 'C': apply H**H (Conjugate transpose)DIRECTDIRECT is CHARACTER*1 Indicates how H is formed from a product of elementary reflectors = 'F': H = H(1) H(2) . . . H(k) (Forward) = 'B': H = H(k) . . . H(2) H(1) (Backward)STOREVSTOREV is CHARACTER*1 Indicates how the vectors which define the elementary reflectors are stored: = 'C': Columns = 'R': RowsMM is INTEGER The number of rows of the matrix B. M >= 0.NN is INTEGER The number of columns of the matrix B. N >= 0.KK is INTEGER The order of the matrix T, i.e. the number of elementary reflectors whose product defines the block reflector. K >= 0.LL is INTEGER The order of the trapezoidal part of V. K >= L >= 0. See Further Details.VV is COMPLEX*16 array, dimension (LDV,K) if STOREV = 'C' (LDV,M) if STOREV = 'R' and SIDE = 'L' (LDV,N) if STOREV = 'R' and SIDE = 'R' The pentagonal matrix V, which contains the elementary reflectors H(1), H(2), ..., H(K). See Further Details.LDVLDV is INTEGER The leading dimension of the array V. If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); if STOREV = 'R', LDV >= K.TT is COMPLEX*16 array, dimension (LDT,K) The triangular K-by-K matrix T in the representation of the block reflector.LDTLDT is INTEGER The leading dimension of the array T. LDT >= K.AA is COMPLEX*16 array, dimension (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' On entry, the K-by-N or M-by-K matrix A. On exit, A is overwritten by the corresponding block of H*C or H**H*C or C*H or C*H**H. See Further Details.LDALDA is INTEGER The leading dimension of the array A. If SIDE = 'L', LDC >= max(1,K); If SIDE = 'R', LDC >= max(1,M).BB is COMPLEX*16 array, dimension (LDB,N) On entry, the M-by-N matrix B. On exit, B is overwritten by the corresponding block of H*C or H**H*C or C*H or C*H**H. See Further Details.LDBLDB is INTEGER The leading dimension of the array B. LDB >= max(1,M).WORKWORK is COMPLEX*16 array, dimension (LDWORK,N) if SIDE = 'L', (LDWORK,K) if SIDE = 'R'.LDWORKLDWORK is INTEGER The leading dimension of the array WORK. If SIDE = 'L', LDWORK >= K; if SIDE = 'R', LDWORK >= M.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:The matrix C is a composite matrix formed from blocks A and B. The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, and if SIDE = 'L', A is of size K-by-N. If SIDE = 'R' and DIRECT = 'F', C = [A B]. If SIDE = 'L' and DIRECT = 'F', C = [A] [B]. If SIDE = 'R' and DIRECT = 'B', C = [B A]. If SIDE = 'L' and DIRECT = 'B', C = [B] [A]. The pentagonal matrix V is composed of a rectangular block V1 and a trapezoidal block V2. The size of the trapezoidal block is determined by the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; if L=0, there is no trapezoidal block, thus V = V1 is rectangular. If DIRECT = 'F' and STOREV = 'C': V = [V1] [V2] - V2 is upper trapezoidal (first L rows of K-by-K upper triangular) If DIRECT = 'F' and STOREV = 'R': V = [V1 V2] - V2 is lower trapezoidal (first L columns of K-by-K lower triangular) If DIRECT = 'B' and STOREV = 'C': V = [V2] [V1] - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K. If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L. If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.

**Author**

Generated automatically by Doxygen for LAPACK from the source code.