Provided by: liblapack-doc_3.7.1-4ubuntu1_all

**NAME**

auxOTHERcomputational

**SYNOPSIS**

Functionscharacter *1 functionchla_transtype(TRANS)CHLA_TRANSTYPEsubroutinedbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)DBDSDCsubroutinedbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)DBDSQRsubroutineddisna(JOB, M, N, D, SEP, INFO)DDISNAsubroutinedlaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)DLAED0used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method. subroutinedlaed1(N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)DLAED1used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal. subroutinedlaed2(K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO)DLAED2used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal. subroutinedlaed3(K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO)DLAED3used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal. subroutinedlaed4(N, I, D, Z, DELTA, RHO, DLAM, INFO)DLAED4used by sstedc. Finds a single root of the secular equation. subroutinedlaed5(I, D, Z, DELTA, RHO, DLAM)DLAED5used by sstedc. Solves the 2-by-2 secular equation. subroutinedlaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)DLAED6used by sstedc. Computes one Newton step in solution of the secular equation. subroutinedlaed7(ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)DLAED7used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense. subroutinedlaed8(ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO)DLAED8used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense. subroutinedlaed9(K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)DLAED9used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense. subroutinedlaeda(N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO)DLAEDAused by sstedc. Computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense. subroutinedlagtf(N, A, LAMBDA, B, C, TOL, D, IN, INFO)DLAGTFcomputes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges. subroutinedlamrg(N1, N2, A, DTRD1, DTRD2, INDEX)DLAMRGcreates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order. subroutinedlartgs(X, Y, SIGMA, CS, SN)DLARTGSgenerates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem. subroutinedlasq1(N, D, E, WORK, INFO)DLASQ1computes the singular values of a real square bidiagonal matrix. Used by sbdsqr. subroutinedlasq2(N, Z, INFO)DLASQ2computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr. subroutinedlasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)DLASQ3checks for deflation, computes a shift and calls dqds. Used by sbdsqr. subroutinedlasq4(I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G)DLASQ4computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. subroutinedlasq5(I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS)DLASQ5computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. subroutinedlasq6(I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2)DLASQ6computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. subroutinedlasrt(ID, N, D, INFO)DLASRTsorts numbers in increasing or decreasing order. subroutinedstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)DSTEBZsubroutinedstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)DSTEDCsubroutinedsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)DSTEQRsubroutinedsterf(N, D, E, INFO)DSTERFinteger functioniladiag(DIAG)ILADIAGinteger functionilaprec(PREC)ILAPRECinteger functionilatrans(TRANS)ILATRANSinteger functionilauplo(UPLO)ILAUPLOsubroutinesbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)SBDSDCsubroutinesbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)SBDSQRsubroutinesdisna(JOB, M, N, D, SEP, INFO)SDISNAsubroutineslaed0(ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, WORK, IWORK, INFO)SLAED0used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method. subroutineslaed1(N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)SLAED1used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal. subroutineslaed2(K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO)SLAED2used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal. subroutineslaed3(K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, CTOT, W, S, INFO)SLAED3used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal. subroutineslaed4(N, I, D, Z, DELTA, RHO, DLAM, INFO)SLAED4used by sstedc. Finds a single root of the secular equation. subroutineslaed5(I, D, Z, DELTA, RHO, DLAM)SLAED5used by sstedc. Solves the 2-by-2 secular equation. subroutineslaed6(KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO)SLAED6used by sstedc. Computes one Newton step in solution of the secular equation. subroutineslaed7(ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)SLAED7used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense. subroutineslaed8(ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, INFO)SLAED8used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense. subroutineslaed9(K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)SLAED9used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense. subroutineslaeda(N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, Q, QPTR, Z, ZTEMP, INFO)SLAEDAused by sstedc. Computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense. subroutineslagtf(N, A, LAMBDA, B, C, TOL, D, IN, INFO)SLAGTFcomputes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges. subroutineslamrg(N1, N2, A, STRD1, STRD2, INDEX)SLAMRGcreates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order. subroutineslartgs(X, Y, SIGMA, CS, SN)SLARTGSgenerates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem. subroutineslasq1(N, D, E, WORK, INFO)SLASQ1computes the singular values of a real square bidiagonal matrix. Used by sbdsqr. subroutineslasq2(N, Z, INFO)SLASQ2computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr. subroutineslasq3(I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, DN2, G, TAU)SLASQ3checks for deflation, computes a shift and calls dqds. Used by sbdsqr. subroutineslasq4(I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, TAU, TTYPE, G)SLASQ4computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. subroutineslasq5(I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, IEEE, EPS)SLASQ5computesonedqdstransforminping-pongform.Usedbysbdsqrandsstegr.subroutineslasq6(I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2)SLASQ6computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. subroutineslasrt(ID, N, D, INFO)SLASRTsorts numbers in increasing or decreasing order. subroutinespttrf(N, D, E, INFO)SPTTRFsubroutinesstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)SSTEBZsubroutinesstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)SSTEDCsubroutinessteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)SSTEQRsubroutinessterf(N, D, E, INFO)SSTERF

**Detailed** **Description**

This is the group of auxiliary Computational routines

**Function** **Documentation**

character*1functionchla_transtype(integerTRANS)CHLA_TRANSTYPEPurpose:This subroutine translates from a BLAST-specified integer constant to the character string specifying a transposition operation. CHLA_TRANSTYPE returns an CHARACTER*1. If CHLA_TRANSTYPE is 'X', then input is not an integer indicating a transposition operator. Otherwise CHLA_TRANSTYPE returns the constant value corresponding to TRANS.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinedbdsdc(characterUPLO,characterCOMPQ,integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,doubleprecision,dimension(ldu,*)U,integerLDU,doubleprecision,dimension(ldvt,*)VT,integerLDVT,doubleprecision,dimension(*)Q,integer,dimension(*)IQ,doubleprecision,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)DBDSDCPurpose:DBDSDC computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, using a divide and conquer method, where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices of left and right singular vectors, respectively. DBDSDC can be used to compute all singular values, and optionally, singular vectors or singular vectors in compact form. This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See DLASD3 for details. The code currently calls DLASDQ if singular values only are desired. However, it can be slightly modified to compute singular values using the divide and conquer method.Parameters:UPLOUPLO is CHARACTER*1 = 'U': B is upper bidiagonal. = 'L': B is lower bidiagonal.COMPQCOMPQ is CHARACTER*1 Specifies whether singular vectors are to be computed as follows: = 'N': Compute singular values only; = 'P': Compute singular values and compute singular vectors in compact form; = 'I': Compute singular values and singular vectors.NN is INTEGER The order of the matrix B. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B.EE is DOUBLE PRECISION array, dimension (N-1) On entry, the elements of E contain the offdiagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E has been destroyed.UU is DOUBLE PRECISION array, dimension (LDU,N) If COMPQ = 'I', then: On exit, if INFO = 0, U contains the left singular vectors of the bidiagonal matrix. For other values of COMPQ, U is not referenced.LDULDU is INTEGER The leading dimension of the array U. LDU >= 1. If singular vectors are desired, then LDU >= max( 1, N ).VTVT is DOUBLE PRECISION array, dimension (LDVT,N) If COMPQ = 'I', then: On exit, if INFO = 0, VT**T contains the right singular vectors of the bidiagonal matrix. For other values of COMPQ, VT is not referenced.LDVTLDVT is INTEGER The leading dimension of the array VT. LDVT >= 1. If singular vectors are desired, then LDVT >= max( 1, N ).QQ is DOUBLE PRECISION array, dimension (LDQ) If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, Q contains all the DOUBLE PRECISION data in LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) words of memory, where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about 25). For other values of COMPQ, Q is not referenced.IQIQ is INTEGER array, dimension (LDIQ) If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, IQ contains all INTEGER data in LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) words of memory, where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about 25). For other values of COMPQ, IQ is not referenced.WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) If COMPQ = 'N' then LWORK >= (4 * N). If COMPQ = 'P' then LWORK >= (6 * N). If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).IWORKIWORK is INTEGER array, dimension (8*N)INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute a singular value. The update process of divide and conquer failed.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USAsubroutinedbdsqr(characterUPLO,integerN,integerNCVT,integerNRU,integerNCC,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,doubleprecision,dimension(ldvt,*)VT,integerLDVT,doubleprecision,dimension(ldu,*)U,integerLDU,doubleprecision,dimension(ldc,*)C,integerLDC,doubleprecision,dimension(*)WORK,integerINFO)DBDSQRPurpose:DBDSQR computes the singular values and, optionally, the right and/or left singular vectors from the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B using the implicit zero-shift QR algorithm. The SVD of B has the form B = Q * S * P**T where S is the diagonal matrix of singular values, Q is an orthogonal matrix of left singular vectors, and P is an orthogonal matrix of right singular vectors. If left singular vectors are requested, this subroutine actually returns U*Q instead of Q, and, if right singular vectors are requested, this subroutine returns P**T*VT instead of P**T, for given real input matrices U and VT. When U and VT are the orthogonal matrices that reduce a general matrix A to bidiagonal form: A = U*B*VT, as computed by DGEBRD, then A = (U*Q) * S * (P**T*VT) is the SVD of A. Optionally, the subroutine may also compute Q**T*C for a given real input matrix C. See "Computing Small Singular Values of Bidiagonal Matrices With Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, no. 5, pp. 873-912, Sept 1990) and "Accurate singular values and differential qd algorithms," by B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics Department, University of California at Berkeley, July 1992 for a detailed description of the algorithm.Parameters:UPLOUPLO is CHARACTER*1 = 'U': B is upper bidiagonal; = 'L': B is lower bidiagonal.NN is INTEGER The order of the matrix B. N >= 0.NCVTNCVT is INTEGER The number of columns of the matrix VT. NCVT >= 0.NRUNRU is INTEGER The number of rows of the matrix U. NRU >= 0.NCCNCC is INTEGER The number of columns of the matrix C. NCC >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B in decreasing order.EE is DOUBLE PRECISION array, dimension (N-1) On entry, the N-1 offdiagonal elements of the bidiagonal matrix B. On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E will contain the diagonal and superdiagonal elements of a bidiagonal matrix orthogonally equivalent to the one given as input.VTVT is DOUBLE PRECISION array, dimension (LDVT, NCVT) On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P**T * VT. Not referenced if NCVT = 0.LDVTLDVT is INTEGER The leading dimension of the array VT. LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.UU is DOUBLE PRECISION array, dimension (LDU, N) On entry, an NRU-by-N matrix U. On exit, U is overwritten by U * Q. Not referenced if NRU = 0.LDULDU is INTEGER The leading dimension of the array U. LDU >= max(1,NRU).CC is DOUBLE PRECISION array, dimension (LDC, NCC) On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q**T * C. Not referenced if NCC = 0.LDCLDC is INTEGER The leading dimension of the array C. LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.WORKWORK is DOUBLE PRECISION array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit < 0: If INFO = -i, the i-th argument had an illegal value > 0: if NCVT = NRU = NCC = 0, = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop) = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks) else NCVT = NRU = NCC = 0, the algorithm did not converge; D and E contain the elements of a bidiagonal matrix which is orthogonally similar to the input matrix B; if INFO = i, i elements of E have not converged to zero.InternalParameters:TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8))) TOLMUL controls the convergence criterion of the QR loop. If it is positive, TOLMUL*EPS is the desired relative precision in the computed singular values. If it is negative, abs(TOLMUL*EPS*sigma_max) is the desired absolute accuracy in the computed singular values (corresponds to relative accuracy abs(TOLMUL*EPS) in the largest singular value. abs(TOLMUL) should be between 1 and 1/EPS, and preferably between 10 (for fast convergence) and .1/EPS (for there to be some accuracy in the results). Default is to lose at either one eighth or 2 of the available decimal digits in each computed singular value (whichever is smaller). MAXITR INTEGER, default = 6 MAXITR controls the maximum number of passes of the algorithm through its inner loop. The algorithms stops (and so fails to converge) if the number of passes through the inner loop exceeds MAXITR*N**2.Note:Bug report from Cezary Dendek. On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is removed since it can overflow pretty easily (for N larger or equal than 18,919). We instead use MAXITDIVN = MAXITR*N.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017subroutineddisna(characterJOB,integerM,integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)SEP,integerINFO)DDISNAPurpose:DDISNA computes the reciprocal condition numbers for the eigenvectors of a real symmetric or complex Hermitian matrix or for the left or right singular vectors of a general m-by-n matrix. The reciprocal condition number is the 'gap' between the corresponding eigenvalue or singular value and the nearest other one. The bound on the error, measured by angle in radians, in the I-th computed vector is given by DLAMCH( 'E' ) * ( ANORM / SEP( I ) ) where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed to be smaller than DLAMCH( 'E' )*ANORM in order to limit the size of the error bound. DDISNA may also be used to compute error bounds for eigenvectors of the generalized symmetric definite eigenproblem.Parameters:JOBJOB is CHARACTER*1 Specifies for which problem the reciprocal condition numbers should be computed: = 'E': the eigenvectors of a symmetric/Hermitian matrix; = 'L': the left singular vectors of a general matrix; = 'R': the right singular vectors of a general matrix.MM is INTEGER The number of rows of the matrix. M >= 0.NN is INTEGER If JOB = 'L' or 'R', the number of columns of the matrix, in which case N >= 0. Ignored if JOB = 'E'.DD is DOUBLE PRECISION array, dimension (M) if JOB = 'E' dimension (min(M,N)) if JOB = 'L' or 'R' The eigenvalues (if JOB = 'E') or singular values (if JOB = 'L' or 'R') of the matrix, in either increasing or decreasing order. If singular values, they must be non-negative.SEPSEP is DOUBLE PRECISION array, dimension (M) if JOB = 'E' dimension (min(M,N)) if JOB = 'L' or 'R' The reciprocal condition numbers of the vectors.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:December 2016subroutinedlaed0(integerICOMPQ,integerQSIZ,integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,doubleprecision,dimension(ldq,*)Q,integerLDQ,doubleprecision,dimension(ldqs,*)QSTORE,integerLDQS,doubleprecision,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)DLAED0used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.Purpose:DLAED0 computes all eigenvalues and corresponding eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method.Parameters:ICOMPQICOMPQ is INTEGER = 0: Compute eigenvalues only. = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form. = 2: Compute eigenvalues and eigenvectors of tridiagonal matrix.QSIZQSIZ is INTEGER The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the main diagonal of the tridiagonal matrix. On exit, its eigenvalues.EE is DOUBLE PRECISION array, dimension (N-1) The off-diagonal elements of the tridiagonal matrix. On exit, E has been destroyed.QQ is DOUBLE PRECISION array, dimension (LDQ, N) On entry, Q must contain an N-by-N orthogonal matrix. If ICOMPQ = 0 Q is not referenced. If ICOMPQ = 1 On entry, Q is a subset of the columns of the orthogonal matrix used to reduce the full matrix to tridiagonal form corresponding to the subset of the full matrix which is being decomposed at this time. If ICOMPQ = 2 On entry, Q will be the identity matrix. On exit, Q contains the eigenvectors of the tridiagonal matrix.LDQLDQ is INTEGER The leading dimension of the array Q. If eigenvectors are desired, then LDQ >= max(1,N). In any case, LDQ >= 1.QSTOREQSTORE is DOUBLE PRECISION array, dimension (LDQS, N) Referenced only when ICOMPQ = 1. Used to store parts of the eigenvector matrix when the updating matrix multiplies take place.LDQSLDQS is INTEGER The leading dimension of the array QSTORE. If ICOMPQ = 1, then LDQS >= max(1,N). In any case, LDQS >= 1.WORKWORK is DOUBLE PRECISION array, If ICOMPQ = 0 or 1, the dimension of WORK must be at least 1 + 3*N + 2*N*lg N + 3*N**2 ( lg( N ) = smallest integer k such that 2^k >= N ) If ICOMPQ = 2, the dimension of WORK must be at least 4*N + N**2.IWORKIWORK is INTEGER array, If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 6 + 6*N + 5*N*lg N. ( lg( N ) = smallest integer k such that 2^k >= N ) If ICOMPQ = 2, the dimension of IWORK must be at least 3 + 5*N.INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutinedlaed1(integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,doubleprecisionRHO,integerCUTPNT,doubleprecision,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)DLAED1used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal.Purpose:DLAED1 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles the case in which eigenvalues only or eigenvalues and eigenvectors of a full symmetric matrix (which was reduced to tridiagonal form) are desired. T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) where Z = Q**T*u, u is a vector of length N with ones in the CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages: The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine DLAED2. The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine DLAED4 (as called by DLAED3). This routine also calculates the eigenvectors of the current problem. The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.Parameters:NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix.QQ is DOUBLE PRECISION array, dimension (LDQ,N) On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) On entry, the permutation which separately sorts the two subproblems in D into ascending order. On exit, the permutation which will reintegrate the subproblems back into sorted order, i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.RHORHO is DOUBLE PRECISION The subdiagonal entry used to create the rank-1 modification.CUTPNTCUTPNT is INTEGER The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N/2.WORKWORK is DOUBLE PRECISION array, dimension (4*N + N**2)IWORKIWORK is INTEGER array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutinedlaed2(integerK,integerN,integerN1,doubleprecision,dimension(*)D,doubleprecision,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,doubleprecisionRHO,doubleprecision,dimension(*)Z,doubleprecision,dimension(*)DLAMDA,doubleprecision,dimension(*)W,doubleprecision,dimension(*)Q2,integer,dimension(*)INDX,integer,dimension(*)INDXC,integer,dimension(*)INDXP,integer,dimension(*)COLTYP,integerINFO)DLAED2used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal.Purpose:DLAED2 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny entry in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one.Parameters:KK is INTEGER The number of non-deflated eigenvalues, and the order of the related secular equation. 0 <= K <=N.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.N1N1 is INTEGER The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= N1 <= N/2.DD is DOUBLE PRECISION array, dimension (N) On entry, D contains the eigenvalues of the two submatrices to be combined. On exit, D contains the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.QQ is DOUBLE PRECISION array, dimension (LDQ, N) On entry, Q contains the eigenvectors of two submatrices in the two square blocks with corners at (1,1), (N1,N1) and (N1+1, N1+1), (N,N). On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have N1 added to their values. Destroyed on exit.RHORHO is DOUBLE PRECISION On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by DLAED3.ZZ is DOUBLE PRECISION array, dimension (N) On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z have been destroyed by the updating process.DLAMDADLAMDA is DOUBLE PRECISION array, dimension (N) A copy of the first K eigenvalues which will be used by DLAED3 to form the secular equation.WW is DOUBLE PRECISION array, dimension (N) The first k values of the final deflation-altered z-vector which will be passed to DLAED3.Q2Q2 is DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2) A copy of the first K eigenvectors which will be used by DLAED3 in a matrix multiply (DGEMM) to solve for the new eigenvectors.INDXINDX is INTEGER array, dimension (N) The permutation used to sort the contents of DLAMDA into ascending order.INDXCINDXC is INTEGER array, dimension (N) The permutation used to arrange the columns of the deflated Q matrix into three groups: the first group contains non-zero elements only at and above N1, the second contains non-zero elements only below N1, and the third is dense.INDXPINDXP is INTEGER array, dimension (N) The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.COLTYPCOLTYP is INTEGER array, dimension (N) During execution, a label which will indicate which of the following types a column in the Q2 matrix is: 1 : non-zero in the upper half only; 2 : dense; 3 : non-zero in the lower half only; 4 : deflated. On exit, COLTYP(i) is the number of columns of type i, for i=1 to 4 only.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:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutinedlaed3(integerK,integerN,integerN1,doubleprecision,dimension(*)D,doubleprecision,dimension(ldq,*)Q,integerLDQ,doubleprecisionRHO,doubleprecision,dimension(*)DLAMDA,doubleprecision,dimension(*)Q2,integer,dimension(*)INDX,integer,dimension(*)CTOT,doubleprecision,dimension(*)W,doubleprecision,dimension(*)S,integerINFO)DLAED3used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal.Purpose:DLAED3 finds the roots of the secular equation, as defined by the values in D, W, and RHO, between 1 and K. It makes the appropriate calls to DLAED4 and then updates the eigenvectors by multiplying the matrix of eigenvectors of the pair of eigensystems being combined by the matrix of eigenvectors of the K-by-K system which is solved here. This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none.Parameters:KK is INTEGER The number of terms in the rational function to be solved by DLAED4. K >= 0.NN is INTEGER The number of rows and columns in the Q matrix. N >= K (deflation may result in N>K).N1N1 is INTEGER The location of the last eigenvalue in the leading submatrix. min(1,N) <= N1 <= N/2.DD is DOUBLE PRECISION array, dimension (N) D(I) contains the updated eigenvalues for 1 <= I <= K.QQ is DOUBLE PRECISION array, dimension (LDQ,N) Initially the first K columns are used as workspace. On output the columns 1 to K contain the updated eigenvectors.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).RHORHO is DOUBLE PRECISION The value of the parameter in the rank one update equation. RHO >= 0 required.DLAMDADLAMDA is DOUBLE PRECISION array, dimension (K) The first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation. May be changed on output by having lowest order bit set to zero on Cray X-MP, Cray Y-MP, Cray-2, or Cray C-90, as described above.Q2Q2 is DOUBLE PRECISION array, dimension (LDQ2*N) The first K columns of this matrix contain the non-deflated eigenvectors for the split problem.INDXINDX is INTEGER array, dimension (N) The permutation used to arrange the columns of the deflated Q matrix into three groups (see DLAED2). The rows of the eigenvectors found by DLAED4 must be likewise permuted before the matrix multiply can take place.CTOTCTOT is INTEGER array, dimension (4) A count of the total number of the various types of columns in Q, as described in INDX. The fourth column type is any column which has been deflated.WW is DOUBLE PRECISION array, dimension (K) The first K elements of this array contain the components of the deflation-adjusted updating vector. Destroyed on output.SS is DOUBLE PRECISION array, dimension (N1 + 1)*K Will contain the eigenvectors of the repaired matrix which will be multiplied by the previously accumulated eigenvectors to update the system.INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutinedlaed4(integerN,integerI,doubleprecision,dimension(*)D,doubleprecision,dimension(*)Z,doubleprecision,dimension(*)DELTA,doubleprecisionRHO,doubleprecisionDLAM,integerINFO)DLAED4used by sstedc. Finds a single root of the secular equation.Purpose:This subroutine computes the I-th updated eigenvalue of a symmetric rank-one modification to a diagonal matrix whose elements are given in the array d, and that D(i) < D(j) for i < j and that RHO > 0. This is arranged by the calling routine, and is no loss in generality. The rank-one modified system is thus diag( D ) + RHO * Z * Z_transpose. where we assume the Euclidean norm of Z is 1. The method consists of approximating the rational functions in the secular equation by simpler interpolating rational functions.Parameters:NN is INTEGER The length of all arrays.II is INTEGER The index of the eigenvalue to be computed. 1 <= I <= N.DD is DOUBLE PRECISION array, dimension (N) The original eigenvalues. It is assumed that they are in order, D(I) < D(J) for I < J.ZZ is DOUBLE PRECISION array, dimension (N) The components of the updating vector.DELTADELTA is DOUBLE PRECISION array, dimension (N) If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5 for detail. The vector DELTA contains the information necessary to construct the eigenvectors by DLAED3 and DLAED9.RHORHO is DOUBLE PRECISION The scalar in the symmetric updating formula.DLAMDLAM is DOUBLE PRECISION The computed lambda_I, the I-th updated eigenvalue.INFOINFO is INTEGER = 0: successful exit > 0: if INFO = 1, the updating process failed.InternalParameters:Logical variable ORGATI (origin-at-i?) is used for distinguishing whether D(i) or D(i+1) is treated as the origin. ORGATI = .true. origin at i ORGATI = .false. origin at i+1 Logical variable SWTCH3 (switch-for-3-poles?) is for noting if we are working with THREE poles! MAXIT is the maximum number of iterations allowed for each eigenvalue.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Ren-Cang Li, Computer Science Division, University of California at Berkeley, USAsubroutinedlaed5(integerI,doubleprecision,dimension(2)D,doubleprecision,dimension(2)Z,doubleprecision,dimension(2)DELTA,doubleprecisionRHO,doubleprecisionDLAM)DLAED5used by sstedc. Solves the 2-by-2 secular equation.Purpose:This subroutine computes the I-th eigenvalue of a symmetric rank-one modification of a 2-by-2 diagonal matrix diag( D ) + RHO * Z * transpose(Z) . The diagonal elements in the array D are assumed to satisfy D(i) < D(j) for i < j . We also assume RHO > 0 and that the Euclidean norm of the vector Z is one.Parameters:II is INTEGER The index of the eigenvalue to be computed. I = 1 or I = 2.DD is DOUBLE PRECISION array, dimension (2) The original eigenvalues. We assume D(1) < D(2).ZZ is DOUBLE PRECISION array, dimension (2) The components of the updating vector.DELTADELTA is DOUBLE PRECISION array, dimension (2) The vector DELTA contains the information necessary to construct the eigenvectors.RHORHO is DOUBLE PRECISION The scalar in the symmetric updating formula.DLAMDLAM is DOUBLE PRECISION The computed lambda_I, the I-th updated eigenvalue.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Ren-Cang Li, Computer Science Division, University of California at Berkeley, USAsubroutinedlaed6(integerKNITER,logicalORGATI,doubleprecisionRHO,doubleprecision,dimension(3)D,doubleprecision,dimension(3)Z,doubleprecisionFINIT,doubleprecisionTAU,integerINFO)DLAED6used by sstedc. Computes one Newton step in solution of the secular equation.Purpose:DLAED6 computes the positive or negative root (closest to the origin) of z(1) z(2) z(3) f(x) = rho + --------- + ---------- + --------- d(1)-x d(2)-x d(3)-x It is assumed that if ORGATI = .true. the root is between d(2) and d(3); otherwise it is between d(1) and d(2) This routine will be called by DLAED4 when necessary. In most cases, the root sought is the smallest in magnitude, though it might not be in some extremely rare situations.Parameters:KNITERKNITER is INTEGER Refer to DLAED4 for its significance.ORGATIORGATI is LOGICAL If ORGATI is true, the needed root is between d(2) and d(3); otherwise it is between d(1) and d(2). See DLAED4 for further details.RHORHO is DOUBLE PRECISION Refer to the equation f(x) above.DD is DOUBLE PRECISION array, dimension (3) D satisfies d(1) < d(2) < d(3).ZZ is DOUBLE PRECISION array, dimension (3) Each of the elements in z must be positive.FINITFINIT is DOUBLE PRECISION The value of f at 0. It is more accurate than the one evaluated inside this routine (if someone wants to do so).TAUTAU is DOUBLE PRECISION The root of the equation f(x).INFOINFO is INTEGER = 0: successful exit > 0: if INFO = 1, failure to convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:10/02/03: This version has a few statements commented out for thread safety (machine parameters are computed on each entry). SJH. 05/10/06: Modified from a new version of Ren-Cang Li, use Gragg-Thornton-Warner cubic convergent scheme for better stability.Contributors:Ren-Cang Li, Computer Science Division, University of California at Berkeley, USAsubroutinedlaed7(integerICOMPQ,integerN,integerQSIZ,integerTLVLS,integerCURLVL,integerCURPBM,doubleprecision,dimension(*)D,doubleprecision,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,doubleprecisionRHO,integerCUTPNT,doubleprecision,dimension(*)QSTORE,integer,dimension(*)QPTR,integer,dimension(*)PRMPTR,integer,dimension(*)PERM,integer,dimension(*)GIVPTR,integer,dimension(2,*)GIVCOL,doubleprecision,dimension(2,*)GIVNUM,doubleprecision,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)DLAED7used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense.Purpose:DLAED7 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and optionally eigenvectors of a dense symmetric matrix that has been reduced to tridiagonal form. DLAED1 handles the case in which all eigenvalues and eigenvectors of a symmetric tridiagonal matrix are desired. T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) where Z = Q**Tu, u is a vector of length N with ones in the CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages: The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine DLAED8. The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine DLAED4 (as called by DLAED9). This routine also calculates the eigenvectors of the current problem. The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.Parameters:ICOMPQICOMPQ is INTEGER = 0: Compute eigenvalues only. = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.QSIZQSIZ is INTEGER The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.TLVLSTLVLS is INTEGER The total number of merging levels in the overall divide and conquer tree.CURLVLCURLVL is INTEGER The current level in the overall merge routine, 0 <= CURLVL <= TLVLS.CURPBMCURPBM is INTEGER The current problem in the current level in the overall merge routine (counting from upper left to lower right).DD is DOUBLE PRECISION array, dimension (N) On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix.QQ is DOUBLE PRECISION array, dimension (LDQ, N) On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) The permutation which will reintegrate the subproblem just solved back into sorted order, i.e., D( INDXQ( I = 1, N ) ) will be in ascending order.RHORHO is DOUBLE PRECISION The subdiagonal element used to create the rank-1 modification.CUTPNTCUTPNT is INTEGER Contains the location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.QSTOREQSTORE is DOUBLE PRECISION array, dimension (N**2+1) Stores eigenvectors of submatrices encountered during divide and conquer, packed together. QPTR points to beginning of the submatrices.QPTRQPTR is INTEGER array, dimension (N+2) List of indices pointing to beginning of submatrices stored in QSTORE. The submatrices are numbered starting at the bottom left of the divide and conquer tree, from left to right and bottom to top.PRMPTRPRMPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and also the size of the full, non-deflated problem.PERMPERM is INTEGER array, dimension (N lg N) Contains the permutations (from deflation and sorting) to be applied to each eigenblock.GIVPTRGIVPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations.GIVCOLGIVCOL is INTEGER array, dimension (2, N lg N) Each pair of numbers indicates a pair of columns to take place in a Givens rotation.GIVNUMGIVNUM is DOUBLE PRECISION array, dimension (2, N lg N) Each number indicates the S value to be used in the corresponding Givens rotation.WORKWORK is DOUBLE PRECISION array, dimension (3*N+2*QSIZ*N)IWORKIWORK is INTEGER array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutinedlaed8(integerICOMPQ,integerK,integerN,integerQSIZ,doubleprecision,dimension(*)D,doubleprecision,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,doubleprecisionRHO,integerCUTPNT,doubleprecision,dimension(*)Z,doubleprecision,dimension(*)DLAMDA,doubleprecision,dimension(ldq2,*)Q2,integerLDQ2,doubleprecision,dimension(*)W,integer,dimension(*)PERM,integerGIVPTR,integer,dimension(2,*)GIVCOL,doubleprecision,dimension(2,*)GIVNUM,integer,dimension(*)INDXP,integer,dimension(*)INDX,integerINFO)DLAED8used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.Purpose:DLAED8 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny element in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one.Parameters:ICOMPQICOMPQ is INTEGER = 0: Compute eigenvalues only. = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.KK is INTEGER The number of non-deflated eigenvalues, and the order of the related secular equation.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.QSIZQSIZ is INTEGER The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.DD is DOUBLE PRECISION array, dimension (N) On entry, the eigenvalues of the two submatrices to be combined. On exit, the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.QQ is DOUBLE PRECISION array, dimension (LDQ,N) If ICOMPQ = 0, Q is not referenced. Otherwise, on entry, Q contains the eigenvectors of the partially solved system which has been previously updated in matrix multiplies with other partially solved eigensystems. On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have CUTPNT added to their values in order to be accurate.RHORHO is DOUBLE PRECISION On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by DLAED3.CUTPNTCUTPNT is INTEGER The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.ZZ is DOUBLE PRECISION array, dimension (N) On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z are destroyed by the updating process.DLAMDADLAMDA is DOUBLE PRECISION array, dimension (N) A copy of the first K eigenvalues which will be used by DLAED3 to form the secular equation.Q2Q2 is DOUBLE PRECISION array, dimension (LDQ2,N) If ICOMPQ = 0, Q2 is not referenced. Otherwise, a copy of the first K eigenvectors which will be used by DLAED7 in a matrix multiply (DGEMM) to update the new eigenvectors.LDQ2LDQ2 is INTEGER The leading dimension of the array Q2. LDQ2 >= max(1,N).WW is DOUBLE PRECISION array, dimension (N) The first k values of the final deflation-altered z-vector and will be passed to DLAED3.PERMPERM is INTEGER array, dimension (N) The permutations (from deflation and sorting) to be applied to each eigenblock.GIVPTRGIVPTR is INTEGER The number of Givens rotations which took place in this subproblem.GIVCOLGIVCOL is INTEGER array, dimension (2, N) Each pair of numbers indicates a pair of columns to take place in a Givens rotation.GIVNUMGIVNUM is DOUBLE PRECISION array, dimension (2, N) Each number indicates the S value to be used in the corresponding Givens rotation.INDXPINDXP is INTEGER array, dimension (N) The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.INDXINDX is INTEGER array, dimension (N) The permutation used to sort the contents of D into ascending order.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:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutinedlaed9(integerK,integerKSTART,integerKSTOP,integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(ldq,*)Q,integerLDQ,doubleprecisionRHO,doubleprecision,dimension(*)DLAMDA,doubleprecision,dimension(*)W,doubleprecision,dimension(lds,*)S,integerLDS,integerINFO)DLAED9used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.Purpose:DLAED9 finds the roots of the secular equation, as defined by the values in D, Z, and RHO, between KSTART and KSTOP. It makes the appropriate calls to DLAED4 and then stores the new matrix of eigenvectors for use in calculating the next level of Z vectors.Parameters:KK is INTEGER The number of terms in the rational function to be solved by DLAED4. K >= 0.KSTARTKSTART is INTEGERKSTOPKSTOP is INTEGER The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP are to be computed. 1 <= KSTART <= KSTOP <= K.NN is INTEGER The number of rows and columns in the Q matrix. N >= K (delation may result in N > K).DD is DOUBLE PRECISION array, dimension (N) D(I) contains the updated eigenvalues for KSTART <= I <= KSTOP.QQ is DOUBLE PRECISION array, dimension (LDQ,N)LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max( 1, N ).RHORHO is DOUBLE PRECISION The value of the parameter in the rank one update equation. RHO >= 0 required.DLAMDADLAMDA is DOUBLE PRECISION array, dimension (K) The first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation.WW is DOUBLE PRECISION array, dimension (K) The first K elements of this array contain the components of the deflation-adjusted updating vector.SS is DOUBLE PRECISION array, dimension (LDS, K) Will contain the eigenvectors of the repaired matrix which will be stored for subsequent Z vector calculation and multiplied by the previously accumulated eigenvectors to update the system.LDSLDS is INTEGER The leading dimension of S. LDS >= max( 1, K ).INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutinedlaeda(integerN,integerTLVLS,integerCURLVL,integerCURPBM,integer,dimension(*)PRMPTR,integer,dimension(*)PERM,integer,dimension(*)GIVPTR,integer,dimension(2,*)GIVCOL,doubleprecision,dimension(2,*)GIVNUM,doubleprecision,dimension(*)Q,integer,dimension(*)QPTR,doubleprecision,dimension(*)Z,doubleprecision,dimension(*)ZTEMP,integerINFO)DLAEDAused by sstedc. Computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense.Purpose:DLAEDA computes the Z vector corresponding to the merge step in the CURLVLth step of the merge process with TLVLS steps for the CURPBMth problem.Parameters:NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.TLVLSTLVLS is INTEGER The total number of merging levels in the overall divide and conquer tree.CURLVLCURLVL is INTEGER The current level in the overall merge routine, 0 <= curlvl <= tlvls.CURPBMCURPBM is INTEGER The current problem in the current level in the overall merge routine (counting from upper left to lower right).PRMPTRPRMPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and incidentally the size of the full, non-deflated problem.PERMPERM is INTEGER array, dimension (N lg N) Contains the permutations (from deflation and sorting) to be applied to each eigenblock.GIVPTRGIVPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations.GIVCOLGIVCOL is INTEGER array, dimension (2, N lg N) Each pair of numbers indicates a pair of columns to take place in a Givens rotation.GIVNUMGIVNUM is DOUBLE PRECISION array, dimension (2, N lg N) Each number indicates the S value to be used in the corresponding Givens rotation.QQ is DOUBLE PRECISION array, dimension (N**2) Contains the square eigenblocks from previous levels, the starting positions for blocks are given by QPTR.QPTRQPTR is INTEGER array, dimension (N+2) Contains a list of pointers which indicate where in Q an eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates the size of the block.ZZ is DOUBLE PRECISION array, dimension (N) On output this vector contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix).ZTEMPZTEMP is DOUBLE PRECISION array, dimension (N)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:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutinedlagtf(integerN,doubleprecision,dimension(*)A,doubleprecisionLAMBDA,doubleprecision,dimension(*)B,doubleprecision,dimension(*)C,doubleprecisionTOL,doubleprecision,dimension(*)D,integer,dimension(*)IN,integerINFO)DLAGTFcomputes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges.Purpose:DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n tridiagonal matrix and lambda is a scalar, as T - lambda*I = PLU, where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column. The factorization is obtained by Gaussian elimination with partial pivoting and implicit row scaling. The parameter LAMBDA is included in the routine so that DLAGTF may be used, in conjunction with DLAGTS, to obtain eigenvectors of T by inverse iteration.Parameters:NN is INTEGER The order of the matrix T.AA is DOUBLE PRECISION array, dimension (N) On entry, A must contain the diagonal elements of T. On exit, A is overwritten by the n diagonal elements of the upper triangular matrix U of the factorization of T.LAMBDALAMBDA is DOUBLE PRECISION On entry, the scalar lambda.BB is DOUBLE PRECISION array, dimension (N-1) On entry, B must contain the (n-1) super-diagonal elements of T. On exit, B is overwritten by the (n-1) super-diagonal elements of the matrix U of the factorization of T.CC is DOUBLE PRECISION array, dimension (N-1) On entry, C must contain the (n-1) sub-diagonal elements of T. On exit, C is overwritten by the (n-1) sub-diagonal elements of the matrix L of the factorization of T.TOLTOL is DOUBLE PRECISION On entry, a relative tolerance used to indicate whether or not the matrix (T - lambda*I) is nearly singular. TOL should normally be chose as approximately the largest relative error in the elements of T. For example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5*10**(-4). If TOL is supplied as less than eps, where eps is the relative machine precision, then the value eps is used in place of TOL.DD is DOUBLE PRECISION array, dimension (N-2) On exit, D is overwritten by the (n-2) second super-diagonal elements of the matrix U of the factorization of T.ININ is INTEGER array, dimension (N) On exit, IN contains details of the permutation matrix P. If an interchange occurred at the kth step of the elimination, then IN(k) = 1, otherwise IN(k) = 0. The element IN(n) returns the smallest positive integer j such that abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL, where norm( A(j) ) denotes the sum of the absolute values of the jth row of the matrix A. If no such j exists then IN(n) is returned as zero. If IN(n) is returned as positive, then a diagonal element of U is small, indicating that (T - lambda*I) is singular or nearly singular,INFOINFO is INTEGER = 0 : successful exit .lt. 0: if INFO = -k, the kth argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinedlamrg(integerN1,integerN2,doubleprecision,dimension(*)A,integerDTRD1,integerDTRD2,integer,dimension(*)INDEX)DLAMRGcreates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order.Purpose:DLAMRG will create a permutation list which will merge the elements of A (which is composed of two independently sorted sets) into a single set which is sorted in ascending order.Parameters:N1N1 is INTEGERN2N2 is INTEGER These arguments contain the respective lengths of the two sorted lists to be merged.AA is DOUBLE PRECISION array, dimension (N1+N2) The first N1 elements of A contain a list of numbers which are sorted in either ascending or descending order. Likewise for the final N2 elements.DTRD1DTRD1 is INTEGERDTRD2DTRD2 is INTEGER These are the strides to be taken through the array A. Allowable strides are 1 and -1. They indicate whether a subset of A is sorted in ascending (DTRDx = 1) or descending (DTRDx = -1) order.INDEXINDEX is INTEGER array, dimension (N1+N2) On exit this array will contain a permutation such that if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be sorted in ascending order.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinedlartgs(doubleprecisionX,doubleprecisionY,doubleprecisionSIGMA,doubleprecisionCS,doubleprecisionSN)DLARTGSgenerates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem.Purpose:DLARTGS generates a plane rotation designed to introduce a bulge in Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD problem. X and Y are the top-row entries, and SIGMA is the shift. The computed CS and SN define a plane rotation satisfying [ CS SN ] . [ X^2 - SIGMA ] = [ R ], [ -SN CS ] [ X * Y ] [ 0 ] with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the rotation is by PI/2.Parameters:XX is DOUBLE PRECISION The (1,1) entry of an upper bidiagonal matrix.YY is DOUBLE PRECISION The (1,2) entry of an upper bidiagonal matrix.SIGMASIGMA is DOUBLE PRECISION The shift.CSCS is DOUBLE PRECISION The cosine of the rotation.SNSN is DOUBLE PRECISION The sine of the rotation.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinedlasq1(integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,doubleprecision,dimension(*)WORK,integerINFO)DLASQ1computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.Purpose:DLASQ1 computes the singular values of a real N-by-N bidiagonal matrix with diagonal D and off-diagonal E. The singular values are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. The algorithm was first presented in "Accurate singular values and differential qd algorithms" by K. V. Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, 1994, and the present implementation is described in "An implementation of the dqds Algorithm (Positive Case)", LAPACK Working Note.Parameters:NN is INTEGER The number of rows and columns in the matrix. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, D contains the diagonal elements of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in decreasing order.EE is DOUBLE PRECISION array, dimension (N) On entry, elements E(1:N-1) contain the off-diagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E is overwritten.WORKWORK is DOUBLE PRECISION array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop) On exit D and E represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into DLASQ1 = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinedlasq2(integerN,doubleprecision,dimension(*)Z,integerINFO)DLASQ2computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.Purpose:DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd array Z to high relative accuracy are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. To see the relation of Z to the tridiagonal matrix, let L be a unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and let U be an upper bidiagonal matrix with 1's above and diagonal Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the symmetric tridiagonal to which it is similar. Note : DLASQ2 defines a logical variable, IEEE, which is true on machines which follow ieee-754 floating-point standard in their handling of infinities and NaNs, and false otherwise. This variable is passed to DLASQ3.Parameters:NN is INTEGER The number of rows and columns in the matrix. N >= 0.ZZ is DOUBLE PRECISION array, dimension ( 4*N ) On entry Z holds the qd array. On exit, entries 1 to N hold the eigenvalues in decreasing order, Z( 2*N+1 ) holds the trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of shifts that failed.INFOINFO is INTEGER = 0: successful exit < 0: if the i-th argument is a scalar and had an illegal value, then INFO = -i, if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j) > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit Z holds a qd array with the same eigenvalues as the given Z. = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:Local Variables: I0:N0 defines a current unreduced segment of Z. The shifts are accumulated in SIGMA. Iteration count is in ITER. Ping-pong is controlled by PP (alternates between 0 and 1).subroutinedlasq3(integerI0,integerN0,doubleprecision,dimension(*)Z,integerPP,doubleprecisionDMIN,doubleprecisionSIGMA,doubleprecisionDESIG,doubleprecisionQMAX,integerNFAIL,integerITER,integerNDIV,logicalIEEE,integerTTYPE,doubleprecisionDMIN1,doubleprecisionDMIN2,doubleprecisionDN,doubleprecisionDN1,doubleprecisionDN2,doubleprecisionG,doubleprecisionTAU)DLASQ3checks for deflation, computes a shift and calls dqds. Used by sbdsqr.Purpose:DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. In case of failure it changes shifts, and tries again until output is positive.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is DOUBLE PRECISION array, dimension ( 4*N0 ) Z holds the qd array.PPPP is INTEGER PP=0 for ping, PP=1 for pong. PP=2 indicates that flipping was applied to the Z array and that the initial tests for deflation should not be performed.DMINDMIN is DOUBLE PRECISION Minimum value of d.SIGMASIGMA is DOUBLE PRECISION Sum of shifts used in current segment.DESIGDESIG is DOUBLE PRECISION Lower order part of SIGMAQMAXQMAX is DOUBLE PRECISION Maximum value of q.NFAILNFAIL is INTEGER Increment NFAIL by 1 each time the shift was too big.ITERITER is INTEGER Increment ITER by 1 for each iteration.NDIVNDIV is INTEGER Increment NDIV by 1 for each division.IEEEIEEE is LOGICAL Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).TTYPETTYPE is INTEGER Shift type.DMIN1DMIN1 is DOUBLE PRECISIONDMIN2DMIN2 is DOUBLE PRECISIONDNDN is DOUBLE PRECISIONDN1DN1 is DOUBLE PRECISIONDN2DN2 is DOUBLE PRECISIONGG is DOUBLE PRECISIONTAUTAU is DOUBLE PRECISION These are passed as arguments in order to save their values between calls to DLASQ3.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinedlasq4(integerI0,integerN0,doubleprecision,dimension(*)Z,integerPP,integerN0IN,doubleprecisionDMIN,doubleprecisionDMIN1,doubleprecisionDMIN2,doubleprecisionDN,doubleprecisionDN1,doubleprecisionDN2,doubleprecisionTAU,integerTTYPE,doubleprecisionG)DLASQ4computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.Purpose:DLASQ4 computes an approximation TAU to the smallest eigenvalue using values of d from the previous transform.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is DOUBLE PRECISION array, dimension ( 4*N0 ) Z holds the qd array.PPPP is INTEGER PP=0 for ping, PP=1 for pong.N0INN0IN is INTEGER The value of N0 at start of EIGTEST.DMINDMIN is DOUBLE PRECISION Minimum value of d.DMIN1DMIN1 is DOUBLE PRECISION Minimum value of d, excluding D( N0 ).DMIN2DMIN2 is DOUBLE PRECISION Minimum value of d, excluding D( N0 ) and D( N0-1 ).DNDN is DOUBLE PRECISION d(N)DN1DN1 is DOUBLE PRECISION d(N-1)DN2DN2 is DOUBLE PRECISION d(N-2)TAUTAU is DOUBLE PRECISION This is the shift.TTYPETTYPE is INTEGER Shift type.GG is DOUBLE PRECISION G is passed as an argument in order to save its value between calls to DLASQ4.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016FurtherDetails:CNST1 = 9/16subroutinedlasq5(integerI0,integerN0,doubleprecision,dimension(*)Z,integerPP,doubleprecisionTAU,doubleprecisionSIGMA,doubleprecisionDMIN,doubleprecisionDMIN1,doubleprecisionDMIN2,doubleprecisionDN,doubleprecisionDNM1,doubleprecisionDNM2,logicalIEEE,doubleprecisionEPS)DLASQ5computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr.Purpose:DLASQ5 computes one dqds transform in ping-pong form, one version for IEEE machines another for non IEEE machines.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is DOUBLE PRECISION array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument.PPPP is INTEGER PP=0 for ping, PP=1 for pong.TAUTAU is DOUBLE PRECISION This is the shift.SIGMASIGMA is DOUBLE PRECISION This is the accumulated shift up to this step.DMINDMIN is DOUBLE PRECISION Minimum value of d.DMIN1DMIN1 is DOUBLE PRECISION Minimum value of d, excluding D( N0 ).DMIN2DMIN2 is DOUBLE PRECISION Minimum value of d, excluding D( N0 ) and D( N0-1 ).DNDN is DOUBLE PRECISION d(N0), the last value of d.DNM1DNM1 is DOUBLE PRECISION d(N0-1).DNM2DNM2 is DOUBLE PRECISION d(N0-2).IEEEIEEE is LOGICAL Flag for IEEE or non IEEE arithmetic.EPSEPS is DOUBLE PRECISION This is the value of epsilon used.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017subroutinedlasq6(integerI0,integerN0,doubleprecision,dimension(*)Z,integerPP,doubleprecisionDMIN,doubleprecisionDMIN1,doubleprecisionDMIN2,doubleprecisionDN,doubleprecisionDNM1,doubleprecisionDNM2)DLASQ6computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.Purpose:DLASQ6 computes one dqd (shift equal to zero) transform in ping-pong form, with protection against underflow and overflow.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is DOUBLE PRECISION array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument.PPPP is INTEGER PP=0 for ping, PP=1 for pong.DMINDMIN is DOUBLE PRECISION Minimum value of d.DMIN1DMIN1 is DOUBLE PRECISION Minimum value of d, excluding D( N0 ).DMIN2DMIN2 is DOUBLE PRECISION Minimum value of d, excluding D( N0 ) and D( N0-1 ).DNDN is DOUBLE PRECISION d(N0), the last value of d.DNM1DNM1 is DOUBLE PRECISION d(N0-1).DNM2DNM2 is DOUBLE PRECISION d(N0-2).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinedlasrt(characterID,integerN,doubleprecision,dimension(*)D,integerINFO)DLASRTsorts numbers in increasing or decreasing order.Purpose:Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ). Use Quick Sort, reverting to Insertion sort on arrays of size <= 20. Dimension of STACK limits N to about 2**32.Parameters:IDID is CHARACTER*1 = 'I': sort D in increasing order; = 'D': sort D in decreasing order.NN is INTEGER The length of the array D.DD is DOUBLE PRECISION array, dimension (N) On entry, the array to be sorted. On exit, D has been sorted into increasing order (D(1) <= ... <= D(N) ) or into decreasing order (D(1) >= ... >= D(N) ), depending on ID.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinedstebz(characterRANGE,characterORDER,integerN,doubleprecisionVL,doubleprecisionVU,integerIL,integerIU,doubleprecisionABSTOL,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,integerM,integerNSPLIT,doubleprecision,dimension(*)W,integer,dimension(*)IBLOCK,integer,dimension(*)ISPLIT,doubleprecision,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)DSTEBZPurpose:DSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix T. The user may ask for all eigenvalues, all eigenvalues in the half-open interval (VL, VU], or the IL-th through IU-th eigenvalues. To avoid overflow, the matrix must be scaled so that its largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest accuracy, it should not be much smaller than that. See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966.Parameters:RANGERANGE is CHARACTER*1 = 'A': ("All") all eigenvalues will be found. = 'V': ("Value") all eigenvalues in the half-open interval (VL, VU] will be found. = 'I': ("Index") the IL-th through IU-th eigenvalues (of the entire matrix) will be found.ORDERORDER is CHARACTER*1 = 'B': ("By Block") the eigenvalues will be grouped by split-off block (see IBLOCK, ISPLIT) and ordered from smallest to largest within the block. = 'E': ("Entire matrix") the eigenvalues for the entire matrix will be ordered from smallest to largest.NN is INTEGER The order of the tridiagonal matrix T. N >= 0.VLVL is DOUBLE PRECISION If RANGE='V', the lower bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.VUVU is DOUBLE PRECISION If RANGE='V', the upper bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.ILIL is INTEGER If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.IUIU is INTEGER If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.ABSTOLABSTOL is DOUBLE PRECISION The absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less. If ABSTOL is less than or equal to zero, then ULP*|T| will be used, where |T| means the 1-norm of T. Eigenvalues will be computed most accurately when ABSTOL is set to twice the underflow threshold 2*DLAMCH('S'), not zero.DD is DOUBLE PRECISION array, dimension (N) The n diagonal elements of the tridiagonal matrix T.EE is DOUBLE PRECISION array, dimension (N-1) The (n-1) off-diagonal elements of the tridiagonal matrix T.MM is INTEGER The actual number of eigenvalues found. 0 <= M <= N. (See also the description of INFO=2,3.)NSPLITNSPLIT is INTEGER The number of diagonal blocks in the matrix T. 1 <= NSPLIT <= N.WW is DOUBLE PRECISION array, dimension (N) On exit, the first M elements of W will contain the eigenvalues. (DSTEBZ may use the remaining N-M elements as workspace.)IBLOCKIBLOCK is INTEGER array, dimension (N) At each row/column j where E(j) is zero or small, the matrix T is considered to split into a block diagonal matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which block (from 1 to the number of blocks) the eigenvalue W(i) belongs. (DSTEBZ may use the remaining N-M elements as workspace.)ISPLITISPLIT is INTEGER array, dimension (N) The splitting points, at which T breaks up into submatrices. The first submatrix consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. (Only the first NSPLIT elements will actually be used, but since the user cannot know a priori what value NSPLIT will have, N words must be reserved for ISPLIT.)WORKWORK is DOUBLE PRECISION array, dimension (4*N)IWORKIWORK is INTEGER array, dimension (3*N)INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: some or all of the eigenvalues failed to converge or were not computed: =1 or 3: Bisection failed to converge for some eigenvalues; these eigenvalues are flagged by a negative block number. The effect is that the eigenvalues may not be as accurate as the absolute and relative tolerances. This is generally caused by unexpectedly inaccurate arithmetic. =2 or 3: RANGE='I' only: Not all of the eigenvalues IL:IU were found. Effect: M < IU+1-IL Cause: non-monotonic arithmetic, causing the Sturm sequence to be non-monotonic. Cure: recalculate, using RANGE='A', and pick out eigenvalues IL:IU. In some cases, increasing the PARAMETER "FUDGE" may make things work. = 4: RANGE='I', and the Gershgorin interval initially used was too small. No eigenvalues were computed. Probable cause: your machine has sloppy floating-point arithmetic. Cure: Increase the PARAMETER "FUDGE", recompile, and try again.InternalParameters:RELFAC DOUBLE PRECISION, default = 2.0e0 The relative tolerance. An interval (a,b] lies within "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), where "ulp" is the machine precision (distance from 1 to the next larger floating point number.) FUDGE DOUBLE PRECISION, default = 2 A "fudge factor" to widen the Gershgorin intervals. Ideally, a value of 1 should work, but on machines with sloppy arithmetic, this needs to be larger. The default for publicly released versions should be large enough to handle the worst machine around. Note that this has no effect on accuracy of the solution.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinedstedc(characterCOMPZ,integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,doubleprecision,dimension(ldz,*)Z,integerLDZ,doubleprecision,dimension(*)WORK,integerLWORK,integer,dimension(*)IWORK,integerLIWORK,integerINFO)DSTEDCPurpose:DSTEDC computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method. The eigenvectors of a full or band real symmetric matrix can also be found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to tridiagonal form. This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See DLAED3 for details.Parameters:COMPZCOMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'I': Compute eigenvectors of tridiagonal matrix also. = 'V': Compute eigenvectors of original dense symmetric matrix also. On entry, Z contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.EE is DOUBLE PRECISION array, dimension (N-1) On entry, the subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.ZZ is DOUBLE PRECISION array, dimension (LDZ,N) On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.LDZLDZ is INTEGER The leading dimension of the array Z. LDZ >= 1. If eigenvectors are desired, then LDZ >= max(1,N).WORKWORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.LWORKLWORK is INTEGER The dimension of the array WORK. If COMPZ = 'N' or N <= 1 then LWORK must be at least 1. If COMPZ = 'V' and N > 1 then LWORK must be at least ( 1 + 3*N + 2*N*lg N + 4*N**2 ), where lg( N ) = smallest integer k such that 2**k >= N. If COMPZ = 'I' and N > 1 then LWORK must be at least ( 1 + 4*N + N**2 ). Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LWORK need only be max(1,2*(N-1)). If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.IWORKIWORK is INTEGER array, dimension (MAX(1,LIWORK)) On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.LIWORKLIWORK is INTEGER The dimension of the array IWORK. If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1. If COMPZ = 'V' and N > 1 then LIWORK must be at least ( 6 + 6*N + 5*N*lg N ). If COMPZ = 'I' and N > 1 then LIWORK must be at least ( 3 + 5*N ). Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LIWORK need only be 1. If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the IWORK array, returns this value as the first entry of the IWORK array, and no error message related to LIWORK is issued by XERBLA.INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutinedsteqr(characterCOMPZ,integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,doubleprecision,dimension(ldz,*)Z,integerLDZ,doubleprecision,dimension(*)WORK,integerINFO)DSTEQRPurpose:DSTEQR computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a full or band symmetric matrix can also be found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to tridiagonal form.Parameters:COMPZCOMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'V': Compute eigenvalues and eigenvectors of the original symmetric matrix. On entry, Z must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form. = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.NN is INTEGER The order of the matrix. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.EE is DOUBLE PRECISION array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.ZZ is DOUBLE PRECISION array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.LDZLDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).WORKWORK is DOUBLE PRECISION array, dimension (max(1,2*N-2)) If COMPZ = 'N', then WORK is not referenced.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is orthogonally similar to the original matrix.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinedsterf(integerN,doubleprecision,dimension(*)D,doubleprecision,dimension(*)E,integerINFO)DSTERFPurpose:DSTERF computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm.Parameters:NN is INTEGER The order of the matrix. N >= 0.DD is DOUBLE PRECISION array, dimension (N) On entry, the n diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.EE is DOUBLE PRECISION array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed to find all of the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016integerfunctioniladiag(characterDIAG)ILADIAGPurpose:This subroutine translated from a character string specifying if a matrix has unit diagonal or not to the relevant BLAST-specified integer constant. ILADIAG returns an INTEGER. If ILADIAG < 0, then the input is not a character indicating a unit or non-unit diagonal. Otherwise ILADIAG returns the constant value corresponding to DIAG.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016integerfunctionilaprec(characterPREC)ILAPRECPurpose:This subroutine translated from a character string specifying an intermediate precision to the relevant BLAST-specified integer constant. ILAPREC returns an INTEGER. If ILAPREC < 0, then the input is not a character indicating a supported intermediate precision. Otherwise ILAPREC returns the constant value corresponding to PREC.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016integerfunctionilatrans(characterTRANS)ILATRANSPurpose:This subroutine translates from a character string specifying a transposition operation to the relevant BLAST-specified integer constant. ILATRANS returns an INTEGER. If ILATRANS < 0, then the input is not a character indicating a transposition operator. Otherwise ILATRANS returns the constant value corresponding to TRANS.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016integerfunctionilauplo(characterUPLO)ILAUPLOPurpose:This subroutine translated from a character string specifying a upper- or lower-triangular matrix to the relevant BLAST-specified integer constant. ILAUPLO returns an INTEGER. If ILAUPLO < 0, then the input is not a character indicating an upper- or lower-triangular matrix. Otherwise ILAUPLO returns the constant value corresponding to UPLO.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinesbdsdc(characterUPLO,characterCOMPQ,integerN,real,dimension(*)D,real,dimension(*)E,real,dimension(ldu,*)U,integerLDU,real,dimension(ldvt,*)VT,integerLDVT,real,dimension(*)Q,integer,dimension(*)IQ,real,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)SBDSDCPurpose:SBDSDC computes the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, using a divide and conquer method, where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices of left and right singular vectors, respectively. SBDSDC can be used to compute all singular values, and optionally, singular vectors or singular vectors in compact form. This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See SLASD3 for details. The code currently calls SLASDQ if singular values only are desired. However, it can be slightly modified to compute singular values using the divide and conquer method.Parameters:UPLOUPLO is CHARACTER*1 = 'U': B is upper bidiagonal. = 'L': B is lower bidiagonal.COMPQCOMPQ is CHARACTER*1 Specifies whether singular vectors are to be computed as follows: = 'N': Compute singular values only; = 'P': Compute singular values and compute singular vectors in compact form; = 'I': Compute singular values and singular vectors.NN is INTEGER The order of the matrix B. N >= 0.DD is REAL array, dimension (N) On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B.EE is REAL array, dimension (N-1) On entry, the elements of E contain the offdiagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E has been destroyed.UU is REAL array, dimension (LDU,N) If COMPQ = 'I', then: On exit, if INFO = 0, U contains the left singular vectors of the bidiagonal matrix. For other values of COMPQ, U is not referenced.LDULDU is INTEGER The leading dimension of the array U. LDU >= 1. If singular vectors are desired, then LDU >= max( 1, N ).VTVT is REAL array, dimension (LDVT,N) If COMPQ = 'I', then: On exit, if INFO = 0, VT**T contains the right singular vectors of the bidiagonal matrix. For other values of COMPQ, VT is not referenced.LDVTLDVT is INTEGER The leading dimension of the array VT. LDVT >= 1. If singular vectors are desired, then LDVT >= max( 1, N ).QQ is REAL array, dimension (LDQ) If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, Q contains all the REAL data in LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) words of memory, where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about 25). For other values of COMPQ, Q is not referenced.IQIQ is INTEGER array, dimension (LDIQ) If COMPQ = 'P', then: On exit, if INFO = 0, Q and IQ contain the left and right singular vectors in a compact form, requiring O(N log N) space instead of 2*N**2. In particular, IQ contains all INTEGER data in LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) words of memory, where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems at the bottom of the computation tree (usually about 25). For other values of COMPQ, IQ is not referenced.WORKWORK is REAL array, dimension (MAX(1,LWORK)) If COMPQ = 'N' then LWORK >= (4 * N). If COMPQ = 'P' then LWORK >= (6 * N). If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).IWORKIWORK is INTEGER array, dimension (8*N)INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute a singular value. The update process of divide and conquer failed.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USAsubroutinesbdsqr(characterUPLO,integerN,integerNCVT,integerNRU,integerNCC,real,dimension(*)D,real,dimension(*)E,real,dimension(ldvt,*)VT,integerLDVT,real,dimension(ldu,*)U,integerLDU,real,dimension(ldc,*)C,integerLDC,real,dimension(*)WORK,integerINFO)SBDSQRPurpose:SBDSQR computes the singular values and, optionally, the right and/or left singular vectors from the singular value decomposition (SVD) of a real N-by-N (upper or lower) bidiagonal matrix B using the implicit zero-shift QR algorithm. The SVD of B has the form B = Q * S * P**T where S is the diagonal matrix of singular values, Q is an orthogonal matrix of left singular vectors, and P is an orthogonal matrix of right singular vectors. If left singular vectors are requested, this subroutine actually returns U*Q instead of Q, and, if right singular vectors are requested, this subroutine returns P**T*VT instead of P**T, for given real input matrices U and VT. When U and VT are the orthogonal matrices that reduce a general matrix A to bidiagonal form: A = U*B*VT, as computed by SGEBRD, then A = (U*Q) * S * (P**T*VT) is the SVD of A. Optionally, the subroutine may also compute Q**T*C for a given real input matrix C. See "Computing Small Singular Values of Bidiagonal Matrices With Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, no. 5, pp. 873-912, Sept 1990) and "Accurate singular values and differential qd algorithms," by B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics Department, University of California at Berkeley, July 1992 for a detailed description of the algorithm.Parameters:UPLOUPLO is CHARACTER*1 = 'U': B is upper bidiagonal; = 'L': B is lower bidiagonal.NN is INTEGER The order of the matrix B. N >= 0.NCVTNCVT is INTEGER The number of columns of the matrix VT. NCVT >= 0.NRUNRU is INTEGER The number of rows of the matrix U. NRU >= 0.NCCNCC is INTEGER The number of columns of the matrix C. NCC >= 0.DD is REAL array, dimension (N) On entry, the n diagonal elements of the bidiagonal matrix B. On exit, if INFO=0, the singular values of B in decreasing order.EE is REAL array, dimension (N-1) On entry, the N-1 offdiagonal elements of the bidiagonal matrix B. On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E will contain the diagonal and superdiagonal elements of a bidiagonal matrix orthogonally equivalent to the one given as input.VTVT is REAL array, dimension (LDVT, NCVT) On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P**T * VT. Not referenced if NCVT = 0.LDVTLDVT is INTEGER The leading dimension of the array VT. LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.UU is REAL array, dimension (LDU, N) On entry, an NRU-by-N matrix U. On exit, U is overwritten by U * Q. Not referenced if NRU = 0.LDULDU is INTEGER The leading dimension of the array U. LDU >= max(1,NRU).CC is REAL array, dimension (LDC, NCC) On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q**T * C. Not referenced if NCC = 0.LDCLDC is INTEGER The leading dimension of the array C. LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.WORKWORK is REAL array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit < 0: If INFO = -i, the i-th argument had an illegal value > 0: if NCVT = NRU = NCC = 0, = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop) = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks) else NCVT = NRU = NCC = 0, the algorithm did not converge; D and E contain the elements of a bidiagonal matrix which is orthogonally similar to the input matrix B; if INFO = i, i elements of E have not converged to zero.InternalParameters:TOLMUL REAL, default = max(10,min(100,EPS**(-1/8))) TOLMUL controls the convergence criterion of the QR loop. If it is positive, TOLMUL*EPS is the desired relative precision in the computed singular values. If it is negative, abs(TOLMUL*EPS*sigma_max) is the desired absolute accuracy in the computed singular values (corresponds to relative accuracy abs(TOLMUL*EPS) in the largest singular value. abs(TOLMUL) should be between 1 and 1/EPS, and preferably between 10 (for fast convergence) and .1/EPS (for there to be some accuracy in the results). Default is to lose at either one eighth or 2 of the available decimal digits in each computed singular value (whichever is smaller). MAXITR INTEGER, default = 6 MAXITR controls the maximum number of passes of the algorithm through its inner loop. The algorithms stops (and so fails to converge) if the number of passes through the inner loop exceeds MAXITR*N**2.Note:Bug report from Cezary Dendek. On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is removed since it can overflow pretty easily (for N larger or equal than 18,919). We instead use MAXITDIVN = MAXITR*N.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017subroutinesdisna(characterJOB,integerM,integerN,real,dimension(*)D,real,dimension(*)SEP,integerINFO)SDISNAPurpose:SDISNA computes the reciprocal condition numbers for the eigenvectors of a real symmetric or complex Hermitian matrix or for the left or right singular vectors of a general m-by-n matrix. The reciprocal condition number is the 'gap' between the corresponding eigenvalue or singular value and the nearest other one. The bound on the error, measured by angle in radians, in the I-th computed vector is given by SLAMCH( 'E' ) * ( ANORM / SEP( I ) ) where ANORM = 2-norm(A) = max( abs( D(j) ) ). SEP(I) is not allowed to be smaller than SLAMCH( 'E' )*ANORM in order to limit the size of the error bound. SDISNA may also be used to compute error bounds for eigenvectors of the generalized symmetric definite eigenproblem.Parameters:JOBJOB is CHARACTER*1 Specifies for which problem the reciprocal condition numbers should be computed: = 'E': the eigenvectors of a symmetric/Hermitian matrix; = 'L': the left singular vectors of a general matrix; = 'R': the right singular vectors of a general matrix.MM is INTEGER The number of rows of the matrix. M >= 0.NN is INTEGER If JOB = 'L' or 'R', the number of columns of the matrix, in which case N >= 0. Ignored if JOB = 'E'.DD is REAL array, dimension (M) if JOB = 'E' dimension (min(M,N)) if JOB = 'L' or 'R' The eigenvalues (if JOB = 'E') or singular values (if JOB = 'L' or 'R') of the matrix, in either increasing or decreasing order. If singular values, they must be non-negative.SEPSEP is REAL array, dimension (M) if JOB = 'E' dimension (min(M,N)) if JOB = 'L' or 'R' The reciprocal condition numbers of the vectors.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:December 2016subroutineslaed0(integerICOMPQ,integerQSIZ,integerN,real,dimension(*)D,real,dimension(*)E,real,dimension(ldq,*)Q,integerLDQ,real,dimension(ldqs,*)QSTORE,integerLDQS,real,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)SLAED0used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.Purpose:SLAED0 computes all eigenvalues and corresponding eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method.Parameters:ICOMPQICOMPQ is INTEGER = 0: Compute eigenvalues only. = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form. = 2: Compute eigenvalues and eigenvectors of tridiagonal matrix.QSIZQSIZ is INTEGER The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.DD is REAL array, dimension (N) On entry, the main diagonal of the tridiagonal matrix. On exit, its eigenvalues.EE is REAL array, dimension (N-1) The off-diagonal elements of the tridiagonal matrix. On exit, E has been destroyed.QQ is REAL array, dimension (LDQ, N) On entry, Q must contain an N-by-N orthogonal matrix. If ICOMPQ = 0 Q is not referenced. If ICOMPQ = 1 On entry, Q is a subset of the columns of the orthogonal matrix used to reduce the full matrix to tridiagonal form corresponding to the subset of the full matrix which is being decomposed at this time. If ICOMPQ = 2 On entry, Q will be the identity matrix. On exit, Q contains the eigenvectors of the tridiagonal matrix.LDQLDQ is INTEGER The leading dimension of the array Q. If eigenvectors are desired, then LDQ >= max(1,N). In any case, LDQ >= 1.QSTOREQSTORE is REAL array, dimension (LDQS, N) Referenced only when ICOMPQ = 1. Used to store parts of the eigenvector matrix when the updating matrix multiplies take place.LDQSLDQS is INTEGER The leading dimension of the array QSTORE. If ICOMPQ = 1, then LDQS >= max(1,N). In any case, LDQS >= 1.WORKWORK is REAL array, If ICOMPQ = 0 or 1, the dimension of WORK must be at least 1 + 3*N + 2*N*lg N + 3*N**2 ( lg( N ) = smallest integer k such that 2^k >= N ) If ICOMPQ = 2, the dimension of WORK must be at least 4*N + N**2.IWORKIWORK is INTEGER array, If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 6 + 6*N + 5*N*lg N. ( lg( N ) = smallest integer k such that 2^k >= N ) If ICOMPQ = 2, the dimension of IWORK must be at least 3 + 5*N.INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutineslaed1(integerN,real,dimension(*)D,real,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,realRHO,integerCUTPNT,real,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)SLAED1used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is tridiagonal.Purpose:SLAED1 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and eigenvectors of a tridiagonal matrix. SLAED7 handles the case in which eigenvalues only or eigenvalues and eigenvectors of a full symmetric matrix (which was reduced to tridiagonal form) are desired. T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) where Z = Q**T*u, u is a vector of length N with ones in the CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages: The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine SLAED2. The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine SLAED4 (as called by SLAED3). This routine also calculates the eigenvectors of the current problem. The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.Parameters:NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.DD is REAL array, dimension (N) On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix.QQ is REAL array, dimension (LDQ,N) On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) On entry, the permutation which separately sorts the two subproblems in D into ascending order. On exit, the permutation which will reintegrate the subproblems back into sorted order, i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.RHORHO is REAL The subdiagonal entry used to create the rank-1 modification.CUTPNTCUTPNT is INTEGER The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N/2.WORKWORK is REAL array, dimension (4*N + N**2)IWORKIWORK is INTEGER array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutineslaed2(integerK,integerN,integerN1,real,dimension(*)D,real,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,realRHO,real,dimension(*)Z,real,dimension(*)DLAMDA,real,dimension(*)W,real,dimension(*)Q2,integer,dimension(*)INDX,integer,dimension(*)INDXC,integer,dimension(*)INDXP,integer,dimension(*)COLTYP,integerINFO)SLAED2used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal.Purpose:SLAED2 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny entry in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one.Parameters:KK is INTEGER The number of non-deflated eigenvalues, and the order of the related secular equation. 0 <= K <=N.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.N1N1 is INTEGER The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= N1 <= N/2.DD is REAL array, dimension (N) On entry, D contains the eigenvalues of the two submatrices to be combined. On exit, D contains the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.QQ is REAL array, dimension (LDQ, N) On entry, Q contains the eigenvectors of two submatrices in the two square blocks with corners at (1,1), (N1,N1) and (N1+1, N1+1), (N,N). On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have N1 added to their values. Destroyed on exit.RHORHO is REAL On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by SLAED3.ZZ is REAL array, dimension (N) On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z have been destroyed by the updating process.DLAMDADLAMDA is REAL array, dimension (N) A copy of the first K eigenvalues which will be used by SLAED3 to form the secular equation.WW is REAL array, dimension (N) The first k values of the final deflation-altered z-vector which will be passed to SLAED3.Q2Q2 is REAL array, dimension (N1**2+(N-N1)**2) A copy of the first K eigenvectors which will be used by SLAED3 in a matrix multiply (SGEMM) to solve for the new eigenvectors.INDXINDX is INTEGER array, dimension (N) The permutation used to sort the contents of DLAMDA into ascending order.INDXCINDXC is INTEGER array, dimension (N) The permutation used to arrange the columns of the deflated Q matrix into three groups: the first group contains non-zero elements only at and above N1, the second contains non-zero elements only below N1, and the third is dense.INDXPINDXP is INTEGER array, dimension (N) The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.COLTYPCOLTYP is INTEGER array, dimension (N) During execution, a label which will indicate which of the following types a column in the Q2 matrix is: 1 : non-zero in the upper half only; 2 : dense; 3 : non-zero in the lower half only; 4 : deflated. On exit, COLTYP(i) is the number of columns of type i, for i=1 to 4 only.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:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutineslaed3(integerK,integerN,integerN1,real,dimension(*)D,real,dimension(ldq,*)Q,integerLDQ,realRHO,real,dimension(*)DLAMDA,real,dimension(*)Q2,integer,dimension(*)INDX,integer,dimension(*)CTOT,real,dimension(*)W,real,dimension(*)S,integerINFO)SLAED3used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal.Purpose:SLAED3 finds the roots of the secular equation, as defined by the values in D, W, and RHO, between 1 and K. It makes the appropriate calls to SLAED4 and then updates the eigenvectors by multiplying the matrix of eigenvectors of the pair of eigensystems being combined by the matrix of eigenvectors of the K-by-K system which is solved here. This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none.Parameters:KK is INTEGER The number of terms in the rational function to be solved by SLAED4. K >= 0.NN is INTEGER The number of rows and columns in the Q matrix. N >= K (deflation may result in N>K).N1N1 is INTEGER The location of the last eigenvalue in the leading submatrix. min(1,N) <= N1 <= N/2.DD is REAL array, dimension (N) D(I) contains the updated eigenvalues for 1 <= I <= K.QQ is REAL array, dimension (LDQ,N) Initially the first K columns are used as workspace. On output the columns 1 to K contain the updated eigenvectors.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).RHORHO is REAL The value of the parameter in the rank one update equation. RHO >= 0 required.DLAMDADLAMDA is REAL array, dimension (K) The first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation. May be changed on output by having lowest order bit set to zero on Cray X-MP, Cray Y-MP, Cray-2, or Cray C-90, as described above.Q2Q2 is REAL array, dimension (LDQ2*N) The first K columns of this matrix contain the non-deflated eigenvectors for the split problem.INDXINDX is INTEGER array, dimension (N) The permutation used to arrange the columns of the deflated Q matrix into three groups (see SLAED2). The rows of the eigenvectors found by SLAED4 must be likewise permuted before the matrix multiply can take place.CTOTCTOT is INTEGER array, dimension (4) A count of the total number of the various types of columns in Q, as described in INDX. The fourth column type is any column which has been deflated.WW is REAL array, dimension (K) The first K elements of this array contain the components of the deflation-adjusted updating vector. Destroyed on output.SS is REAL array, dimension (N1 + 1)*K Will contain the eigenvectors of the repaired matrix which will be multiplied by the previously accumulated eigenvectors to update the system.INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2017Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutineslaed4(integerN,integerI,real,dimension(*)D,real,dimension(*)Z,real,dimension(*)DELTA,realRHO,realDLAM,integerINFO)SLAED4used by sstedc. Finds a single root of the secular equation.Purpose:This subroutine computes the I-th updated eigenvalue of a symmetric rank-one modification to a diagonal matrix whose elements are given in the array d, and that D(i) < D(j) for i < j and that RHO > 0. This is arranged by the calling routine, and is no loss in generality. The rank-one modified system is thus diag( D ) + RHO * Z * Z_transpose. where we assume the Euclidean norm of Z is 1. The method consists of approximating the rational functions in the secular equation by simpler interpolating rational functions.Parameters:NN is INTEGER The length of all arrays.II is INTEGER The index of the eigenvalue to be computed. 1 <= I <= N.DD is REAL array, dimension (N) The original eigenvalues. It is assumed that they are in order, D(I) < D(J) for I < J.ZZ is REAL array, dimension (N) The components of the updating vector.DELTADELTA is REAL array, dimension (N) If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5 for detail. The vector DELTA contains the information necessary to construct the eigenvectors by SLAED3 and SLAED9.RHORHO is REAL The scalar in the symmetric updating formula.DLAMDLAM is REAL The computed lambda_I, the I-th updated eigenvalue.INFOINFO is INTEGER = 0: successful exit > 0: if INFO = 1, the updating process failed.InternalParameters:Logical variable ORGATI (origin-at-i?) is used for distinguishing whether D(i) or D(i+1) is treated as the origin. ORGATI = .true. origin at i ORGATI = .false. origin at i+1 Logical variable SWTCH3 (switch-for-3-poles?) is for noting if we are working with THREE poles! MAXIT is the maximum number of iterations allowed for each eigenvalue.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Ren-Cang Li, Computer Science Division, University of California at Berkeley, USAsubroutineslaed5(integerI,real,dimension(2)D,real,dimension(2)Z,real,dimension(2)DELTA,realRHO,realDLAM)SLAED5used by sstedc. Solves the 2-by-2 secular equation.Purpose:This subroutine computes the I-th eigenvalue of a symmetric rank-one modification of a 2-by-2 diagonal matrix diag( D ) + RHO * Z * transpose(Z) . The diagonal elements in the array D are assumed to satisfy D(i) < D(j) for i < j . We also assume RHO > 0 and that the Euclidean norm of the vector Z is one.Parameters:II is INTEGER The index of the eigenvalue to be computed. I = 1 or I = 2.DD is REAL array, dimension (2) The original eigenvalues. We assume D(1) < D(2).ZZ is REAL array, dimension (2) The components of the updating vector.DELTADELTA is REAL array, dimension (2) The vector DELTA contains the information necessary to construct the eigenvectors.RHORHO is REAL The scalar in the symmetric updating formula.DLAMDLAM is REAL The computed lambda_I, the I-th updated eigenvalue.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Ren-Cang Li, Computer Science Division, University of California at Berkeley, USAsubroutineslaed6(integerKNITER,logicalORGATI,realRHO,real,dimension(3)D,real,dimension(3)Z,realFINIT,realTAU,integerINFO)SLAED6used by sstedc. Computes one Newton step in solution of the secular equation.Purpose:SLAED6 computes the positive or negative root (closest to the origin) of z(1) z(2) z(3) f(x) = rho + --------- + ---------- + --------- d(1)-x d(2)-x d(3)-x It is assumed that if ORGATI = .true. the root is between d(2) and d(3); otherwise it is between d(1) and d(2) This routine will be called by SLAED4 when necessary. In most cases, the root sought is the smallest in magnitude, though it might not be in some extremely rare situations.Parameters:KNITERKNITER is INTEGER Refer to SLAED4 for its significance.ORGATIORGATI is LOGICAL If ORGATI is true, the needed root is between d(2) and d(3); otherwise it is between d(1) and d(2). See SLAED4 for further details.RHORHO is REAL Refer to the equation f(x) above.DD is REAL array, dimension (3) D satisfies d(1) < d(2) < d(3).ZZ is REAL array, dimension (3) Each of the elements in z must be positive.FINITFINIT is REAL The value of f at 0. It is more accurate than the one evaluated inside this routine (if someone wants to do so).TAUTAU is REAL The root of the equation f(x).INFOINFO is INTEGER = 0: successful exit > 0: if INFO = 1, failure to convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:10/02/03: This version has a few statements commented out for thread safety (machine parameters are computed on each entry). SJH. 05/10/06: Modified from a new version of Ren-Cang Li, use Gragg-Thornton-Warner cubic convergent scheme for better stability.Contributors:Ren-Cang Li, Computer Science Division, University of California at Berkeley, USAsubroutineslaed7(integerICOMPQ,integerN,integerQSIZ,integerTLVLS,integerCURLVL,integerCURPBM,real,dimension(*)D,real,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,realRHO,integerCUTPNT,real,dimension(*)QSTORE,integer,dimension(*)QPTR,integer,dimension(*)PRMPTR,integer,dimension(*)PERM,integer,dimension(*)GIVPTR,integer,dimension(2,*)GIVCOL,real,dimension(2,*)GIVNUM,real,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)SLAED7used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense.Purpose:SLAED7 computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. This routine is used only for the eigenproblem which requires all eigenvalues and optionally eigenvectors of a dense symmetric matrix that has been reduced to tridiagonal form. SLAED1 handles the case in which all eigenvalues and eigenvectors of a symmetric tridiagonal matrix are desired. T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) where Z = Q**Tu, u is a vector of length N with ones in the CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. The eigenvectors of the original matrix are stored in Q, and the eigenvalues are in D. The algorithm consists of three stages: The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in the Z vector. For each such occurrence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine SLAED8. The second stage consists of calculating the updated eigenvalues. This is done by finding the roots of the secular equation via the routine SLAED4 (as called by SLAED9). This routine also calculates the eigenvectors of the current problem. The final stage consists of computing the updated eigenvectors directly using the updated eigenvalues. The eigenvectors for the current problem are multiplied with the eigenvectors from the overall problem.Parameters:ICOMPQICOMPQ is INTEGER = 0: Compute eigenvalues only. = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.QSIZQSIZ is INTEGER The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.TLVLSTLVLS is INTEGER The total number of merging levels in the overall divide and conquer tree.CURLVLCURLVL is INTEGER The current level in the overall merge routine, 0 <= CURLVL <= TLVLS.CURPBMCURPBM is INTEGER The current problem in the current level in the overall merge routine (counting from upper left to lower right).DD is REAL array, dimension (N) On entry, the eigenvalues of the rank-1-perturbed matrix. On exit, the eigenvalues of the repaired matrix.QQ is REAL array, dimension (LDQ, N) On entry, the eigenvectors of the rank-1-perturbed matrix. On exit, the eigenvectors of the repaired tridiagonal matrix.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) The permutation which will reintegrate the subproblem just solved back into sorted order, i.e., D( INDXQ( I = 1, N ) ) will be in ascending order.RHORHO is REAL The subdiagonal element used to create the rank-1 modification.CUTPNTCUTPNT is INTEGER Contains the location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.QSTOREQSTORE is REAL array, dimension (N**2+1) Stores eigenvectors of submatrices encountered during divide and conquer, packed together. QPTR points to beginning of the submatrices.QPTRQPTR is INTEGER array, dimension (N+2) List of indices pointing to beginning of submatrices stored in QSTORE. The submatrices are numbered starting at the bottom left of the divide and conquer tree, from left to right and bottom to top.PRMPTRPRMPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and also the size of the full, non-deflated problem.PERMPERM is INTEGER array, dimension (N lg N) Contains the permutations (from deflation and sorting) to be applied to each eigenblock.GIVPTRGIVPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations.GIVCOLGIVCOL is INTEGER array, dimension (2, N lg N) Each pair of numbers indicates a pair of columns to take place in a Givens rotation.GIVNUMGIVNUM is REAL array, dimension (2, N lg N) Each number indicates the S value to be used in the corresponding Givens rotation.WORKWORK is REAL array, dimension (3*N+2*QSIZ*N)IWORKIWORK is INTEGER array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutineslaed8(integerICOMPQ,integerK,integerN,integerQSIZ,real,dimension(*)D,real,dimension(ldq,*)Q,integerLDQ,integer,dimension(*)INDXQ,realRHO,integerCUTPNT,real,dimension(*)Z,real,dimension(*)DLAMDA,real,dimension(ldq2,*)Q2,integerLDQ2,real,dimension(*)W,integer,dimension(*)PERM,integerGIVPTR,integer,dimension(2,*)GIVCOL,real,dimension(2,*)GIVNUM,integer,dimension(*)INDXP,integer,dimension(*)INDX,integerINFO)SLAED8used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.Purpose:SLAED8 merges the two sets of eigenvalues together into a single sorted set. Then it tries to deflate the size of the problem. There are two ways in which deflation can occur: when two or more eigenvalues are close together or if there is a tiny element in the Z vector. For each such occurrence the order of the related secular equation problem is reduced by one.Parameters:ICOMPQICOMPQ is INTEGER = 0: Compute eigenvalues only. = 1: Compute eigenvectors of original dense symmetric matrix also. On entry, Q contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.KK is INTEGER The number of non-deflated eigenvalues, and the order of the related secular equation.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.QSIZQSIZ is INTEGER The dimension of the orthogonal matrix used to reduce the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.DD is REAL array, dimension (N) On entry, the eigenvalues of the two submatrices to be combined. On exit, the trailing (N-K) updated eigenvalues (those which were deflated) sorted into increasing order.QQ is REAL array, dimension (LDQ,N) If ICOMPQ = 0, Q is not referenced. Otherwise, on entry, Q contains the eigenvectors of the partially solved system which has been previously updated in matrix multiplies with other partially solved eigensystems. On exit, Q contains the trailing (N-K) updated eigenvectors (those which were deflated) in its last N-K columns.LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N).INDXQINDXQ is INTEGER array, dimension (N) The permutation which separately sorts the two sub-problems in D into ascending order. Note that elements in the second half of this permutation must first have CUTPNT added to their values in order to be accurate.RHORHO is REAL On entry, the off-diagonal element associated with the rank-1 cut which originally split the two submatrices which are now being recombined. On exit, RHO has been modified to the value required by SLAED3.CUTPNTCUTPNT is INTEGER The location of the last eigenvalue in the leading sub-matrix. min(1,N) <= CUTPNT <= N.ZZ is REAL array, dimension (N) On entry, Z contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix). On exit, the contents of Z are destroyed by the updating process.DLAMDADLAMDA is REAL array, dimension (N) A copy of the first K eigenvalues which will be used by SLAED3 to form the secular equation.Q2Q2 is REAL array, dimension (LDQ2,N) If ICOMPQ = 0, Q2 is not referenced. Otherwise, a copy of the first K eigenvectors which will be used by SLAED7 in a matrix multiply (SGEMM) to update the new eigenvectors.LDQ2LDQ2 is INTEGER The leading dimension of the array Q2. LDQ2 >= max(1,N).WW is REAL array, dimension (N) The first k values of the final deflation-altered z-vector and will be passed to SLAED3.PERMPERM is INTEGER array, dimension (N) The permutations (from deflation and sorting) to be applied to each eigenblock.GIVPTRGIVPTR is INTEGER The number of Givens rotations which took place in this subproblem.GIVCOLGIVCOL is INTEGER array, dimension (2, N) Each pair of numbers indicates a pair of columns to take place in a Givens rotation.GIVNUMGIVNUM is REAL array, dimension (2, N) Each number indicates the S value to be used in the corresponding Givens rotation.INDXPINDXP is INTEGER array, dimension (N) The permutation used to place deflated values of D at the end of the array. INDXP(1:K) points to the nondeflated D-values and INDXP(K+1:N) points to the deflated eigenvalues.INDXINDX is INTEGER array, dimension (N) The permutation used to sort the contents of D into ascending order.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:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutineslaed9(integerK,integerKSTART,integerKSTOP,integerN,real,dimension(*)D,real,dimension(ldq,*)Q,integerLDQ,realRHO,real,dimension(*)DLAMDA,real,dimension(*)W,real,dimension(lds,*)S,integerLDS,integerINFO)SLAED9used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.Purpose:SLAED9 finds the roots of the secular equation, as defined by the values in D, Z, and RHO, between KSTART and KSTOP. It makes the appropriate calls to SLAED4 and then stores the new matrix of eigenvectors for use in calculating the next level of Z vectors.Parameters:KK is INTEGER The number of terms in the rational function to be solved by SLAED4. K >= 0.KSTARTKSTART is INTEGERKSTOPKSTOP is INTEGER The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP are to be computed. 1 <= KSTART <= KSTOP <= K.NN is INTEGER The number of rows and columns in the Q matrix. N >= K (delation may result in N > K).DD is REAL array, dimension (N) D(I) contains the updated eigenvalues for KSTART <= I <= KSTOP.QQ is REAL array, dimension (LDQ,N)LDQLDQ is INTEGER The leading dimension of the array Q. LDQ >= max( 1, N ).RHORHO is REAL The value of the parameter in the rank one update equation. RHO >= 0 required.DLAMDADLAMDA is REAL array, dimension (K) The first K elements of this array contain the old roots of the deflated updating problem. These are the poles of the secular equation.WW is REAL array, dimension (K) The first K elements of this array contain the components of the deflation-adjusted updating vector.SS is REAL array, dimension (LDS, K) Will contain the eigenvectors of the repaired matrix which will be stored for subsequent Z vector calculation and multiplied by the previously accumulated eigenvectors to update the system.LDSLDS is INTEGER The leading dimension of S. LDS >= max( 1, K ).INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, an eigenvalue did not convergeAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutineslaeda(integerN,integerTLVLS,integerCURLVL,integerCURPBM,integer,dimension(*)PRMPTR,integer,dimension(*)PERM,integer,dimension(*)GIVPTR,integer,dimension(2,*)GIVCOL,real,dimension(2,*)GIVNUM,real,dimension(*)Q,integer,dimension(*)QPTR,real,dimension(*)Z,real,dimension(*)ZTEMP,integerINFO)SLAEDAused by sstedc. Computes the Z vector determining the rank-one modification of the diagonal matrix. Used when the original matrix is dense.Purpose:SLAEDA computes the Z vector corresponding to the merge step in the CURLVLth step of the merge process with TLVLS steps for the CURPBMth problem.Parameters:NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.TLVLSTLVLS is INTEGER The total number of merging levels in the overall divide and conquer tree.CURLVLCURLVL is INTEGER The current level in the overall merge routine, 0 <= curlvl <= tlvls.CURPBMCURPBM is INTEGER The current problem in the current level in the overall merge routine (counting from upper left to lower right).PRMPTRPRMPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in PERM a level's permutation is stored. PRMPTR(i+1) - PRMPTR(i) indicates the size of the permutation and incidentally the size of the full, non-deflated problem.PERMPERM is INTEGER array, dimension (N lg N) Contains the permutations (from deflation and sorting) to be applied to each eigenblock.GIVPTRGIVPTR is INTEGER array, dimension (N lg N) Contains a list of pointers which indicate where in GIVCOL a level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i) indicates the number of Givens rotations.GIVCOLGIVCOL is INTEGER array, dimension (2, N lg N) Each pair of numbers indicates a pair of columns to take place in a Givens rotation.GIVNUMGIVNUM is REAL array, dimension (2, N lg N) Each number indicates the S value to be used in the corresponding Givens rotation.QQ is REAL array, dimension (N**2) Contains the square eigenblocks from previous levels, the starting positions for blocks are given by QPTR.QPTRQPTR is INTEGER array, dimension (N+2) Contains a list of pointers which indicate where in Q an eigenblock is stored. SQRT( QPTR(i+1) - QPTR(i) ) indicates the size of the block.ZZ is REAL array, dimension (N) On output this vector contains the updating vector (the last row of the first sub-eigenvector matrix and the first row of the second sub-eigenvector matrix).ZTEMPZTEMP is REAL array, dimension (N)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:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USAsubroutineslagtf(integerN,real,dimension(*)A,realLAMBDA,real,dimension(*)B,real,dimension(*)C,realTOL,real,dimension(*)D,integer,dimension(*)IN,integerINFO)SLAGTFcomputes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix, and λ a scalar, using partial pivoting with row interchanges.Purpose:SLAGTF factorizes the matrix (T - lambda*I), where T is an n by n tridiagonal matrix and lambda is a scalar, as T - lambda*I = PLU, where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column. The factorization is obtained by Gaussian elimination with partial pivoting and implicit row scaling. The parameter LAMBDA is included in the routine so that SLAGTF may be used, in conjunction with SLAGTS, to obtain eigenvectors of T by inverse iteration.Parameters:NN is INTEGER The order of the matrix T.AA is REAL array, dimension (N) On entry, A must contain the diagonal elements of T. On exit, A is overwritten by the n diagonal elements of the upper triangular matrix U of the factorization of T.LAMBDALAMBDA is REAL On entry, the scalar lambda.BB is REAL array, dimension (N-1) On entry, B must contain the (n-1) super-diagonal elements of T. On exit, B is overwritten by the (n-1) super-diagonal elements of the matrix U of the factorization of T.CC is REAL array, dimension (N-1) On entry, C must contain the (n-1) sub-diagonal elements of T. On exit, C is overwritten by the (n-1) sub-diagonal elements of the matrix L of the factorization of T.TOLTOL is REAL On entry, a relative tolerance used to indicate whether or not the matrix (T - lambda*I) is nearly singular. TOL should normally be chose as approximately the largest relative error in the elements of T. For example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5*10**(-4). If TOL is supplied as less than eps, where eps is the relative machine precision, then the value eps is used in place of TOL.DD is REAL array, dimension (N-2) On exit, D is overwritten by the (n-2) second super-diagonal elements of the matrix U of the factorization of T.ININ is INTEGER array, dimension (N) On exit, IN contains details of the permutation matrix P. If an interchange occurred at the kth step of the elimination, then IN(k) = 1, otherwise IN(k) = 0. The element IN(n) returns the smallest positive integer j such that abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL, where norm( A(j) ) denotes the sum of the absolute values of the jth row of the matrix A. If no such j exists then IN(n) is returned as zero. If IN(n) is returned as positive, then a diagonal element of U is small, indicating that (T - lambda*I) is singular or nearly singular,INFOINFO is INTEGER = 0 : successful exit .lt. 0: if INFO = -k, the kth argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutineslamrg(integerN1,integerN2,real,dimension(*)A,integerSTRD1,integerSTRD2,integer,dimension(*)INDEX)SLAMRGcreates a permutation list to merge the entries of two independently sorted sets into a single set sorted in ascending order.Purpose:SLAMRG will create a permutation list which will merge the elements of A (which is composed of two independently sorted sets) into a single set which is sorted in ascending order.Parameters:N1N1 is INTEGERN2N2 is INTEGER These arguments contain the respective lengths of the two sorted lists to be merged.AA is REAL array, dimension (N1+N2) The first N1 elements of A contain a list of numbers which are sorted in either ascending or descending order. Likewise for the final N2 elements.STRD1STRD1 is INTEGERSTRD2STRD2 is INTEGER These are the strides to be taken through the array A. Allowable strides are 1 and -1. They indicate whether a subset of A is sorted in ascending (STRDx = 1) or descending (STRDx = -1) order.INDEXINDEX is INTEGER array, dimension (N1+N2) On exit this array will contain a permutation such that if B( I ) = A( INDEX( I ) ) for I=1,N1+N2, then B will be sorted in ascending order.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutineslartgs(realX,realY,realSIGMA,realCS,realSN)SLARTGSgenerates a plane rotation designed to introduce a bulge in implicit QR iteration for the bidiagonal SVD problem.Purpose:SLARTGS generates a plane rotation designed to introduce a bulge in Golub-Reinsch-style implicit QR iteration for the bidiagonal SVD problem. X and Y are the top-row entries, and SIGMA is the shift. The computed CS and SN define a plane rotation satisfying [ CS SN ] . [ X^2 - SIGMA ] = [ R ], [ -SN CS ] [ X * Y ] [ 0 ] with R nonnegative. If X^2 - SIGMA and X * Y are 0, then the rotation is by PI/2.Parameters:XX is REAL The (1,1) entry of an upper bidiagonal matrix.YY is REAL The (1,2) entry of an upper bidiagonal matrix.SIGMASIGMA is REAL The shift.CSCS is REAL The cosine of the rotation.SNSN is REAL The sine of the rotation.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutineslasq1(integerN,real,dimension(*)D,real,dimension(*)E,real,dimension(*)WORK,integerINFO)SLASQ1computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.Purpose:SLASQ1 computes the singular values of a real N-by-N bidiagonal matrix with diagonal D and off-diagonal E. The singular values are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. The algorithm was first presented in "Accurate singular values and differential qd algorithms" by K. V. Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, 1994, and the present implementation is described in "An implementation of the dqds Algorithm (Positive Case)", LAPACK Working Note.Parameters:NN is INTEGER The number of rows and columns in the matrix. N >= 0.DD is REAL array, dimension (N) On entry, D contains the diagonal elements of the bidiagonal matrix whose SVD is desired. On normal exit, D contains the singular values in decreasing order.EE is REAL array, dimension (N) On entry, elements E(1:N-1) contain the off-diagonal elements of the bidiagonal matrix whose SVD is desired. On exit, E is overwritten.WORKWORK is REAL array, dimension (4*N)INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop) On exit D and E represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into SLASQ1 = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutineslasq2(integerN,real,dimension(*)Z,integerINFO)SLASQ2computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.Purpose:SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd array Z to high relative accuracy are computed to high relative accuracy, in the absence of denormalization, underflow and overflow. To see the relation of Z to the tridiagonal matrix, let L be a unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and let U be an upper bidiagonal matrix with 1's above and diagonal Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the symmetric tridiagonal to which it is similar. Note : SLASQ2 defines a logical variable, IEEE, which is true on machines which follow ieee-754 floating-point standard in their handling of infinities and NaNs, and false otherwise. This variable is passed to SLASQ3.Parameters:NN is INTEGER The number of rows and columns in the matrix. N >= 0.ZZ is REAL array, dimension ( 4*N ) On entry Z holds the qd array. On exit, entries 1 to N hold the eigenvalues in decreasing order, Z( 2*N+1 ) holds the trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of shifts that failed.INFOINFO is INTEGER = 0: successful exit < 0: if the i-th argument is a scalar and had an illegal value, then INFO = -i, if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j) > 0: the algorithm failed = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 100*N iterations (in inner while loop). On exit Z holds a qd array with the same eigenvalues as the given Z. = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks)Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016FurtherDetails:Local Variables: I0:N0 defines a current unreduced segment of Z. The shifts are accumulated in SIGMA. Iteration count is in ITER. Ping-pong is controlled by PP (alternates between 0 and 1).subroutineslasq3(integerI0,integerN0,real,dimension(*)Z,integerPP,realDMIN,realSIGMA,realDESIG,realQMAX,integerNFAIL,integerITER,integerNDIV,logicalIEEE,integerTTYPE,realDMIN1,realDMIN2,realDN,realDN1,realDN2,realG,realTAU)SLASQ3checks for deflation, computes a shift and calls dqds. Used by sbdsqr.Purpose:SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. In case of failure it changes shifts, and tries again until output is positive.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is REAL array, dimension ( 4*N0 ) Z holds the qd array.PPPP is INTEGER PP=0 for ping, PP=1 for pong. PP=2 indicates that flipping was applied to the Z array and that the initial tests for deflation should not be performed.DMINDMIN is REAL Minimum value of d.SIGMASIGMA is REAL Sum of shifts used in current segment.DESIGDESIG is REAL Lower order part of SIGMAQMAXQMAX is REAL Maximum value of q.NFAILNFAIL is INTEGER Increment NFAIL by 1 each time the shift was too big.ITERITER is INTEGER Increment ITER by 1 for each iteration.NDIVNDIV is INTEGER Increment NDIV by 1 for each division.IEEEIEEE is LOGICAL Flag for IEEE or non IEEE arithmetic (passed to SLASQ5).TTYPETTYPE is INTEGER Shift type.DMIN1DMIN1 is REALDMIN2DMIN2 is REALDNDN is REALDN1DN1 is REALDN2DN2 is REALGG is REALTAUTAU is REAL These are passed as arguments in order to save their values between calls to SLASQ3.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutineslasq4(integerI0,integerN0,real,dimension(*)Z,integerPP,integerN0IN,realDMIN,realDMIN1,realDMIN2,realDN,realDN1,realDN2,realTAU,integerTTYPE,realG)SLASQ4computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr.Purpose:SLASQ4 computes an approximation TAU to the smallest eigenvalue using values of d from the previous transform.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is REAL array, dimension ( 4*N0 ) Z holds the qd array.PPPP is INTEGER PP=0 for ping, PP=1 for pong.N0INN0IN is INTEGER The value of N0 at start of EIGTEST.DMINDMIN is REAL Minimum value of d.DMIN1DMIN1 is REAL Minimum value of d, excluding D( N0 ).DMIN2DMIN2 is REAL Minimum value of d, excluding D( N0 ) and D( N0-1 ).DNDN is REAL d(N)DN1DN1 is REAL d(N-1)DN2DN2 is REAL d(N-2)TAUTAU is REAL This is the shift.TTYPETTYPE is INTEGER Shift type.GG is REAL G is passed as an argument in order to save its value between calls to SLASQ4.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016FurtherDetails:CNST1 = 9/16subroutineslasq5(integerI0,integerN0,real,dimension(*)Z,integerPP,realTAU,realSIGMA,realDMIN,realDMIN1,realDMIN2,realDN,realDNM1,realDNM2,logicalIEEE,realEPS)SLASQ5computesonedqdstransforminping-pongform.Usedbysbdsqrandsstegr.Purpose:SLASQ5 computes one dqds transform in ping-pong form, one version for IEEE machines another for non IEEE machines.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is REAL array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument.PPPP is INTEGER PP=0 for ping, PP=1 for pong.TAUTAU is REAL This is the shift.SIGMASIGMA is REAL This is the accumulated shift up to this step.DMINDMIN is REAL Minimum value of d.DMIN1DMIN1 is REAL Minimum value of d, excluding D( N0 ).DMIN2DMIN2 is REAL Minimum value of d, excluding D( N0 ) and D( N0-1 ).DNDN is REAL d(N0), the last value of d.DNM1DNM1 is REAL d(N0-1).DNM2DNM2 is REAL d(N0-2).IEEEIEEE is LOGICAL Flag for IEEE or non IEEE arithmetic.EPSEPS is REAL This is the value of epsilon used.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutineslasq6(integerI0,integerN0,real,dimension(*)Z,integerPP,realDMIN,realDMIN1,realDMIN2,realDN,realDNM1,realDNM2)SLASQ6computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.Purpose:SLASQ6 computes one dqd (shift equal to zero) transform in ping-pong form, with protection against underflow and overflow.Parameters:I0I0 is INTEGER First index.N0N0 is INTEGER Last index.ZZ is REAL array, dimension ( 4*N ) Z holds the qd array. EMIN is stored in Z(4*N0) to avoid an extra argument.PPPP is INTEGER PP=0 for ping, PP=1 for pong.DMINDMIN is REAL Minimum value of d.DMIN1DMIN1 is REAL Minimum value of d, excluding D( N0 ).DMIN2DMIN2 is REAL Minimum value of d, excluding D( N0 ) and D( N0-1 ).DNDN is REAL d(N0), the last value of d.DNM1DNM1 is REAL d(N0-1).DNM2DNM2 is REAL d(N0-2).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutineslasrt(characterID,integerN,real,dimension(*)D,integerINFO)SLASRTsorts numbers in increasing or decreasing order.Purpose:Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ). Use Quick Sort, reverting to Insertion sort on arrays of size <= 20. Dimension of STACK limits N to about 2**32.Parameters:IDID is CHARACTER*1 = 'I': sort D in increasing order; = 'D': sort D in decreasing order.NN is INTEGER The length of the array D.DD is REAL array, dimension (N) On entry, the array to be sorted. On exit, D has been sorted into increasing order (D(1) <= ... <= D(N) ) or into decreasing order (D(1) >= ... >= D(N) ), depending on ID.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal valueAuthor:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinespttrf(integerN,real,dimension(*)D,real,dimension(*)E,integerINFO)SPTTRFPurpose:SPTTRF computes the L*D*L**T factorization of a real symmetric positive definite tridiagonal matrix A. The factorization may also be regarded as having the form A = U**T*D*U.Parameters:NN is INTEGER The order of the matrix A. N >= 0.DD is REAL array, dimension (N) On entry, the n diagonal elements of the tridiagonal matrix A. On exit, the n diagonal elements of the diagonal matrix D from the L*D*L**T factorization of A.EE is REAL array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix A. On exit, the (n-1) subdiagonal elements of the unit bidiagonal factor L from the L*D*L**T factorization of A. E can also be regarded as the superdiagonal of the unit bidiagonal factor U from the U**T*D*U factorization of A.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal value > 0: if INFO = k, the leading minor of order k is not positive definite; if k < N, the factorization could not be completed, while if k = N, the factorization was completed, but D(N) <= 0.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinesstebz(characterRANGE,characterORDER,integerN,realVL,realVU,integerIL,integerIU,realABSTOL,real,dimension(*)D,real,dimension(*)E,integerM,integerNSPLIT,real,dimension(*)W,integer,dimension(*)IBLOCK,integer,dimension(*)ISPLIT,real,dimension(*)WORK,integer,dimension(*)IWORK,integerINFO)SSTEBZPurpose:SSTEBZ computes the eigenvalues of a symmetric tridiagonal matrix T. The user may ask for all eigenvalues, all eigenvalues in the half-open interval (VL, VU], or the IL-th through IU-th eigenvalues. To avoid overflow, the matrix must be scaled so that its largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest accuracy, it should not be much smaller than that. See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix", Report CS41, Computer Science Dept., Stanford University, July 21, 1966.Parameters:RANGERANGE is CHARACTER*1 = 'A': ("All") all eigenvalues will be found. = 'V': ("Value") all eigenvalues in the half-open interval (VL, VU] will be found. = 'I': ("Index") the IL-th through IU-th eigenvalues (of the entire matrix) will be found.ORDERORDER is CHARACTER*1 = 'B': ("By Block") the eigenvalues will be grouped by split-off block (see IBLOCK, ISPLIT) and ordered from smallest to largest within the block. = 'E': ("Entire matrix") the eigenvalues for the entire matrix will be ordered from smallest to largest.NN is INTEGER The order of the tridiagonal matrix T. N >= 0.VLVL is REAL If RANGE='V', the lower bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.VUVU is REAL If RANGE='V', the upper bound of the interval to be searched for eigenvalues. Eigenvalues less than or equal to VL, or greater than VU, will not be returned. VL < VU. Not referenced if RANGE = 'A' or 'I'.ILIL is INTEGER If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.IUIU is INTEGER If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. Not referenced if RANGE = 'A' or 'V'.ABSTOLABSTOL is REAL The absolute tolerance for the eigenvalues. An eigenvalue (or cluster) is considered to be located if it has been determined to lie in an interval whose width is ABSTOL or less. If ABSTOL is less than or equal to zero, then ULP*|T| will be used, where |T| means the 1-norm of T. Eigenvalues will be computed most accurately when ABSTOL is set to twice the underflow threshold 2*SLAMCH('S'), not zero.DD is REAL array, dimension (N) The n diagonal elements of the tridiagonal matrix T.EE is REAL array, dimension (N-1) The (n-1) off-diagonal elements of the tridiagonal matrix T.MM is INTEGER The actual number of eigenvalues found. 0 <= M <= N. (See also the description of INFO=2,3.)NSPLITNSPLIT is INTEGER The number of diagonal blocks in the matrix T. 1 <= NSPLIT <= N.WW is REAL array, dimension (N) On exit, the first M elements of W will contain the eigenvalues. (SSTEBZ may use the remaining N-M elements as workspace.)IBLOCKIBLOCK is INTEGER array, dimension (N) At each row/column j where E(j) is zero or small, the matrix T is considered to split into a block diagonal matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which block (from 1 to the number of blocks) the eigenvalue W(i) belongs. (SSTEBZ may use the remaining N-M elements as workspace.)ISPLITISPLIT is INTEGER array, dimension (N) The splitting points, at which T breaks up into submatrices. The first submatrix consists of rows/columns 1 to ISPLIT(1), the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), etc., and the NSPLIT-th consists of rows/columns ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. (Only the first NSPLIT elements will actually be used, but since the user cannot know a priori what value NSPLIT will have, N words must be reserved for ISPLIT.)WORKWORK is REAL array, dimension (4*N)IWORKIWORK is INTEGER array, dimension (3*N)INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: some or all of the eigenvalues failed to converge or were not computed: =1 or 3: Bisection failed to converge for some eigenvalues; these eigenvalues are flagged by a negative block number. The effect is that the eigenvalues may not be as accurate as the absolute and relative tolerances. This is generally caused by unexpectedly inaccurate arithmetic. =2 or 3: RANGE='I' only: Not all of the eigenvalues IL:IU were found. Effect: M < IU+1-IL Cause: non-monotonic arithmetic, causing the Sturm sequence to be non-monotonic. Cure: recalculate, using RANGE='A', and pick out eigenvalues IL:IU. In some cases, increasing the PARAMETER "FUDGE" may make things work. = 4: RANGE='I', and the Gershgorin interval initially used was too small. No eigenvalues were computed. Probable cause: your machine has sloppy floating-point arithmetic. Cure: Increase the PARAMETER "FUDGE", recompile, and try again.InternalParameters:RELFAC REAL, default = 2.0e0 The relative tolerance. An interval (a,b] lies within "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), where "ulp" is the machine precision (distance from 1 to the next larger floating point number.) FUDGE REAL, default = 2 A "fudge factor" to widen the Gershgorin intervals. Ideally, a value of 1 should work, but on machines with sloppy arithmetic, this needs to be larger. The default for publicly released versions should be large enough to handle the worst machine around. Note that this has no effect on accuracy of the solution.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:June 2016subroutinesstedc(characterCOMPZ,integerN,real,dimension(*)D,real,dimension(*)E,real,dimension(ldz,*)Z,integerLDZ,real,dimension(*)WORK,integerLWORK,integer,dimension(*)IWORK,integerLIWORK,integerINFO)SSTEDCPurpose:SSTEDC computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the divide and conquer method. The eigenvectors of a full or band real symmetric matrix can also be found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to tridiagonal form. This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract, or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none. See SLAED3 for details.Parameters:COMPZCOMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'I': Compute eigenvectors of tridiagonal matrix also. = 'V': Compute eigenvectors of original dense symmetric matrix also. On entry, Z contains the orthogonal matrix used to reduce the original matrix to tridiagonal form.NN is INTEGER The dimension of the symmetric tridiagonal matrix. N >= 0.DD is REAL array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.EE is REAL array, dimension (N-1) On entry, the subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.ZZ is REAL array, dimension (LDZ,N) On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.LDZLDZ is INTEGER The leading dimension of the array Z. LDZ >= 1. If eigenvectors are desired, then LDZ >= max(1,N).WORKWORK is REAL array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.LWORKLWORK is INTEGER The dimension of the array WORK. If COMPZ = 'N' or N <= 1 then LWORK must be at least 1. If COMPZ = 'V' and N > 1 then LWORK must be at least ( 1 + 3*N + 2*N*lg N + 4*N**2 ), where lg( N ) = smallest integer k such that 2**k >= N. If COMPZ = 'I' and N > 1 then LWORK must be at least ( 1 + 4*N + N**2 ). Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LWORK need only be max(1,2*(N-1)). If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.IWORKIWORK is INTEGER array, dimension (MAX(1,LIWORK)) On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.LIWORKLIWORK is INTEGER The dimension of the array IWORK. If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1. If COMPZ = 'V' and N > 1 then LIWORK must be at least ( 6 + 6*N + 5*N*lg N ). If COMPZ = 'I' and N > 1 then LIWORK must be at least ( 3 + 5*N ). Note that for COMPZ = 'I' or 'V', then if N is less than or equal to the minimum divide size, usually 25, then LIWORK need only be 1. If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the IWORK array, returns this value as the first entry of the IWORK array, and no error message related to LIWORK is issued by XERBLA.INFOINFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: The algorithm failed to compute an eigenvalue while working on the submatrix lying in rows and columns INFO/(N+1) through mod(INFO,N+1).Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016Contributors:Jeff Rutter, Computer Science Division, University of California at Berkeley, USA Modified by Francoise Tisseur, University of Tennesseesubroutinessteqr(characterCOMPZ,integerN,real,dimension(*)D,real,dimension(*)E,real,dimension(ldz,*)Z,integerLDZ,real,dimension(*)WORK,integerINFO)SSTEQRPurpose:SSTEQR computes all eigenvalues and, optionally, eigenvectors of a symmetric tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a full or band symmetric matrix can also be found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to tridiagonal form.Parameters:COMPZCOMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'V': Compute eigenvalues and eigenvectors of the original symmetric matrix. On entry, Z must contain the orthogonal matrix used to reduce the original matrix to tridiagonal form. = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.NN is INTEGER The order of the matrix. N >= 0.DD is REAL array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.EE is REAL array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.ZZ is REAL array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the orthogonal matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original symmetric matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.LDZLDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).WORKWORK is REAL array, dimension (max(1,2*N-2)) If COMPZ = 'N', then WORK is not referenced.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is orthogonally similar to the original matrix.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016subroutinessterf(integerN,real,dimension(*)D,real,dimension(*)E,integerINFO)SSTERFPurpose:SSTERF computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm.Parameters:NN is INTEGER The order of the matrix. N >= 0.DD is REAL array, dimension (N) On entry, the n diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.EE is REAL array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.INFOINFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm failed to find all of the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero.Author:Univ. of Tennessee Univ. of California Berkeley Univ. of Colorado Denver NAG Ltd.Date:December 2016

**Author**

Generated automatically by Doxygen for LAPACK from the source code.