CLATRS3 solves one of the triangular systems


    A * X = B * diag(scale),  A**T * X = B * diag(scale), or


    A**H * X = B * diag(scale)


 with scaling to prevent overflow.  Here A is an upper or lower


 triangular matrix, A**T denotes the transpose of A, A**H denotes the


 conjugate transpose of A. X and B are n-by-nrhs matrices and scale


 is an nrhs-element vector of scaling factors. A scaling factor scale(j)


 is usually less than or equal to 1, chosen such that X(:,j) is less


 than the overflow threshold. If the matrix A is singular (A(j,j) = 0


 for some j), then a non-trivial solution to A*X = 0 is returned. If


 the system is so badly scaled that the solution cannot be represented


 as (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned.


 This is a BLAS-3 version of LATRS for solving several right


 hand sides simultaneously.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


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

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

NORMIN

          NORMIN is CHARACTER*1


          Specifies whether CNORM has been set or not.


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


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


                  be computed and stored in CNORM.

N

          N is INTEGER


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

NRHS

          NRHS is INTEGER


          The number of columns of X.  NRHS >= 0.

A

          A is COMPLEX array, dimension (LDA,N)


          The triangular matrix A.  If UPLO = 'U', the leading n by n


          upper triangular part of the array A contains the upper


          triangular matrix, and the strictly lower triangular part of


          A is not referenced.  If UPLO = 'L', the leading n by n lower


          triangular part of the array A contains the lower triangular


          matrix, and the strictly upper triangular part of A is not


          referenced.  If DIAG = 'U', the diagonal elements of A are


          also not referenced and are assumed to be 1.

LDA

          LDA is INTEGER


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

X

          X is COMPLEX array, dimension (LDX,NRHS)


          On entry, the right hand side B of the triangular system.


          On exit, X is overwritten by the solution matrix X.

LDX

          LDX is INTEGER


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

SCALE

          SCALE is REAL array, dimension (NRHS)


          The scaling factor s(k) is for the triangular system


          A * x(:,k) = s(k)*b(:,k)  or  A**T* x(:,k) = s(k)*b(:,k).


          If SCALE = 0, the matrix A is singular or badly scaled.


          If A(j,j) = 0 is encountered, a non-trivial vector x(:,k)


          that is an exact or approximate solution to A*x(:,k) = 0


          is returned. If the system so badly scaled that solution


          cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0


          is returned.

CNORM

          CNORM is REAL array, dimension (N)


          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)


          contains the norm of the off-diagonal part of the j-th column


          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal


          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)


          must be greater than or equal to the 1-norm.


          If NORMIN = 'N', CNORM is an output argument and CNORM(j)


          returns the 1-norm of the offdiagonal part of the j-th column


          of A.

WORK

          WORK is REAL array, dimension (LWORK).


          On exit, if INFO = 0, WORK(1) returns the optimal size of


          WORK.

LWORK LWORK is INTEGER LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLABRD reduces the first NB rows and columns of a real general


 m by n matrix A to upper or lower bidiagonal form by an orthogonal


 transformation Q**T * A * P, and returns the matrices X and Y which


 are needed to apply the transformation to the unreduced part of A.


 If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower


 bidiagonal form.


 This is an auxiliary routine called by DGEBRD

M

          M is INTEGER


          The number of rows in the matrix A.

N

          N is INTEGER


          The number of columns in the matrix A.

NB

          NB is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the m by n general matrix to be reduced.


          On exit, the first NB rows and columns of the matrix are


          overwritten; the rest of the array is unchanged.


          If m >= n, elements on and below the diagonal in the first NB


            columns, with the array TAUQ, represent the orthogonal


            matrix Q as a product of elementary reflectors; and


            elements above the diagonal in the first NB rows, with the


            array TAUP, represent the orthogonal matrix P as a product


            of elementary reflectors.


          If m < n, elements below the diagonal in the first NB


            columns, with the array TAUQ, represent the orthogonal


            matrix Q as a product of elementary reflectors, and


            elements on and above the diagonal in the first NB rows,


            with the array TAUP, represent the orthogonal matrix P as


            a product of elementary reflectors.


          See Further Details.

LDA

          LDA is INTEGER


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

D

          D is DOUBLE PRECISION array, dimension (NB)


          The diagonal elements of the first NB rows and columns of


          the reduced matrix.  D(i) = A(i,i).

E

          E is DOUBLE PRECISION array, dimension (NB)


          The off-diagonal elements of the first NB rows and columns of


          the reduced matrix.

TAUQ

          TAUQ is DOUBLE PRECISION array, dimension (NB)


          The scalar factors of the elementary reflectors which


          represent the orthogonal matrix Q. See Further Details.

TAUP

          TAUP is DOUBLE PRECISION array, dimension (NB)


          The scalar factors of the elementary reflectors which


          represent the orthogonal matrix P. See Further Details.

X

          X is DOUBLE PRECISION array, dimension (LDX,NB)


          The m-by-nb matrix X required to update the unreduced part


          of A.

LDX

          LDX is INTEGER


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

Y

          Y is DOUBLE PRECISION array, dimension (LDY,NB)


          The n-by-nb matrix Y required to update the unreduced part


          of A.

LDY

          LDY is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The matrices Q and P are represented as products of elementary


  reflectors:


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


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


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


  where tauq and taup are real scalars, and v and u are real vectors.


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


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


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


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


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


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


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


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


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


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


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


  with nb = 2:


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


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


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


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


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


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


    (  v1  v2  a   a   a  )


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


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


  of the vector defining G(i).

 DLACN2 estimates the 1-norm of a square, real matrix A.


 Reverse communication is used for evaluating matrix-vector products.

N

          N is INTEGER


         The order of the matrix.  N >= 1.

V

          V is DOUBLE PRECISION array, dimension (N)


         On the final return, V = A*W,  where  EST = norm(V)/norm(W)


         (W is not returned).

X

          X is DOUBLE PRECISION array, dimension (N)


         On an intermediate return, X should be overwritten by


               A * X,   if KASE=1,


               A**T * X,  if KASE=2,


         and DLACN2 must be re-called with all the other parameters


         unchanged.

ISGN

          ISGN is INTEGER array, dimension (N)

EST

          EST is DOUBLE PRECISION


         On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be


         unchanged from the previous call to DLACN2.


         On exit, EST is an estimate (a lower bound) for norm(A).

KASE

          KASE is INTEGER


         On the initial call to DLACN2, KASE should be 0.


         On an intermediate return, KASE will be 1 or 2, indicating


         whether X should be overwritten by A * X  or A**T * X.


         On the final return from DLACN2, KASE will again be 0.

ISAVE

          ISAVE is INTEGER array, dimension (3)


         ISAVE is used to save variables between calls to DLACN2

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Originally named SONEST, dated March 16, 1988.


  This is a thread safe version of DLACON, which uses the array ISAVE


  in place of a SAVE statement, as follows:


     DLACON     DLACN2


      JUMP     ISAVE(1)


      J        ISAVE(2)


      ITER     ISAVE(3)

Nick Higham, University of Manchester

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

 DLACON estimates the 1-norm of a square, real matrix A.


 Reverse communication is used for evaluating matrix-vector products.

N

          N is INTEGER


         The order of the matrix.  N >= 1.

V

          V is DOUBLE PRECISION array, dimension (N)


         On the final return, V = A*W,  where  EST = norm(V)/norm(W)


         (W is not returned).

X

          X is DOUBLE PRECISION array, dimension (N)


         On an intermediate return, X should be overwritten by


               A * X,   if KASE=1,


               A**T * X,  if KASE=2,


         and DLACON must be re-called with all the other parameters


         unchanged.

ISGN

          ISGN is INTEGER array, dimension (N)

EST

          EST is DOUBLE PRECISION


         On entry with KASE = 1 or 2 and JUMP = 3, EST should be


         unchanged from the previous call to DLACON.


         On exit, EST is an estimate (a lower bound) for norm(A).

KASE

          KASE is INTEGER


         On the initial call to DLACON, KASE should be 0.


         On an intermediate return, KASE will be 1 or 2, indicating


         whether X should be overwritten by A * X  or A**T * X.


         On the final return from DLACON, KASE will again be 0.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Nick Higham, University of Manchester.
Originally named SONEST, dated March 16, 1988.

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

 DLADIV performs complex division in  real arithmetic


                       a + i*b


            p + i*q = ---------


                       c + i*d


 The algorithm is due to Michael Baudin and Robert L. Smith


 and can be found in the paper


 'A Robust Complex Division in Scilab'

A

          A is DOUBLE PRECISION

B

          B is DOUBLE PRECISION

C

          C is DOUBLE PRECISION

D

          D is DOUBLE PRECISION


          The scalars a, b, c, and d in the above expression.

P

          P is DOUBLE PRECISION

Q

          Q is DOUBLE PRECISION


          The scalars p and q in the above expression.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAEIN uses inverse iteration to find a right or left eigenvector


 corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg


 matrix H.

RIGHTV

          RIGHTV is LOGICAL


          = .TRUE. : compute right eigenvector;


          = .FALSE.: compute left eigenvector.

NOINIT

          NOINIT is LOGICAL


          = .TRUE. : no initial vector supplied in (VR,VI).


          = .FALSE.: initial vector supplied in (VR,VI).

N

          N is INTEGER


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

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


          The upper Hessenberg matrix H.

LDH

          LDH is INTEGER


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

WR

          WR is DOUBLE PRECISION

WI

          WI is DOUBLE PRECISION


          The real and imaginary parts of the eigenvalue of H whose


          corresponding right or left eigenvector is to be computed.

VR

          VR is DOUBLE PRECISION array, dimension (N)

VI

          VI is DOUBLE PRECISION array, dimension (N)


          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain


          a real starting vector for inverse iteration using the real


          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI


          must contain the real and imaginary parts of a complex


          starting vector for inverse iteration using the complex


          eigenvalue (WR,WI); otherwise VR and VI need not be set.


          On exit, if WI = 0.0 (real eigenvalue), VR contains the


          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),


          VR and VI contain the real and imaginary parts of the


          computed complex eigenvector. The eigenvector is normalized


          so that the component of largest magnitude has magnitude 1;


          here the magnitude of a complex number (x,y) is taken to be


          |x| + |y|.


          VI is not referenced if WI = 0.0.

B

          B is DOUBLE PRECISION array, dimension (LDB,N)

LDB

          LDB is INTEGER


          The leading dimension of the array B.  LDB >= N+1.

WORK

          WORK is DOUBLE PRECISION array, dimension (N)

EPS3

          EPS3 is DOUBLE PRECISION


          A small machine-dependent value which is used to perturb


          close eigenvalues, and to replace zero pivots.

SMLNUM

          SMLNUM is DOUBLE PRECISION


          A machine-dependent value close to the underflow threshold.

BIGNUM

          BIGNUM is DOUBLE PRECISION


          A machine-dependent value close to the overflow threshold.

INFO

          INFO is INTEGER


          = 0:  successful exit


          = 1:  inverse iteration did not converge; VR is set to the


                last iterate, and so is VI if WI.ne.0.0.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in


 an upper quasi-triangular matrix T by an orthogonal similarity


 transformation.


 T must be in Schur canonical form, that is, block upper triangular


 with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block


 has its diagonal elements equal and its off-diagonal elements of


 opposite sign.

WANTQ

          WANTQ is LOGICAL


          = .TRUE. : accumulate the transformation in the matrix Q;


          = .FALSE.: do not accumulate the transformation.

N

          N is INTEGER


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

T

          T is DOUBLE PRECISION array, dimension (LDT,N)


          On entry, the upper quasi-triangular matrix T, in Schur


          canonical form.


          On exit, the updated matrix T, again in Schur canonical form.

LDT

          LDT is INTEGER


          The leading dimension of the array T. LDT >= max(1,N).

Q

          Q is DOUBLE PRECISION array, dimension (LDQ,N)


          On entry, if WANTQ is .TRUE., the orthogonal matrix Q.


          On exit, if WANTQ is .TRUE., the updated matrix Q.


          If WANTQ is .FALSE., Q is not referenced.

LDQ

          LDQ is INTEGER


          The leading dimension of the array Q.


          LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.

J1

          J1 is INTEGER


          The index of the first row of the first block T11.

N1

          N1 is INTEGER


          The order of the first block T11. N1 = 0, 1 or 2.

N2

          N2 is INTEGER


          The order of the second block T22. N2 = 0, 1 or 2.

WORK

          WORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          = 0: successful exit


          = 1: the transformed matrix T would be too far from Schur


               form; the blocks are not swapped and T and Q are


               unchanged.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue


 problem  A - w B, with scaling as necessary to avoid over-/underflow.


 The scaling factor 's' results in a modified eigenvalue equation


     s A - w B


 where  s  is a non-negative scaling factor chosen so that  w,  w B,


 and  s A  do not overflow and, if possible, do not underflow, either.

A

          A is DOUBLE PRECISION array, dimension (LDA, 2)


          On entry, the 2 x 2 matrix A.  It is assumed that its 1-norm


          is less than 1/SAFMIN.  Entries less than


          sqrt(SAFMIN)*norm(A) are subject to being treated as zero.

LDA

          LDA is INTEGER


          The leading dimension of the array A.  LDA >= 2.

B

          B is DOUBLE PRECISION array, dimension (LDB, 2)


          On entry, the 2 x 2 upper triangular matrix B.  It is


          assumed that the one-norm of B is less than 1/SAFMIN.  The


          diagonals should be at least sqrt(SAFMIN) times the largest


          element of B (in absolute value); if a diagonal is smaller


          than that, then  +/- sqrt(SAFMIN) will be used instead of


          that diagonal.

LDB

          LDB is INTEGER


          The leading dimension of the array B.  LDB >= 2.

SAFMIN

          SAFMIN is DOUBLE PRECISION


          The smallest positive number s.t. 1/SAFMIN does not


          overflow.  (This should always be DLAMCH('S') -- it is an


          argument in order to avoid having to call DLAMCH frequently.)

SCALE1

          SCALE1 is DOUBLE PRECISION


          A scaling factor used to avoid over-/underflow in the


          eigenvalue equation which defines the first eigenvalue.  If


          the eigenvalues are complex, then the eigenvalues are


          ( WR1  +/-  WI i ) / SCALE1  (which may lie outside the


          exponent range of the machine), SCALE1=SCALE2, and SCALE1


          will always be positive.  If the eigenvalues are real, then


          the first (real) eigenvalue is  WR1 / SCALE1 , but this may


          overflow or underflow, and in fact, SCALE1 may be zero or


          less than the underflow threshold if the exact eigenvalue


          is sufficiently large.

SCALE2

          SCALE2 is DOUBLE PRECISION


          A scaling factor used to avoid over-/underflow in the


          eigenvalue equation which defines the second eigenvalue.  If


          the eigenvalues are complex, then SCALE2=SCALE1.  If the


          eigenvalues are real, then the second (real) eigenvalue is


          WR2 / SCALE2 , but this may overflow or underflow, and in


          fact, SCALE2 may be zero or less than the underflow


          threshold if the exact eigenvalue is sufficiently large.

WR1

          WR1 is DOUBLE PRECISION


          If the eigenvalue is real, then WR1 is SCALE1 times the


          eigenvalue closest to the (2,2) element of A B**(-1).  If the


          eigenvalue is complex, then WR1=WR2 is SCALE1 times the real


          part of the eigenvalues.

WR2

          WR2 is DOUBLE PRECISION


          If the eigenvalue is real, then WR2 is SCALE2 times the


          other eigenvalue.  If the eigenvalue is complex, then


          WR1=WR2 is SCALE1 times the real part of the eigenvalues.

WI

          WI is DOUBLE PRECISION


          If the eigenvalue is real, then WI is zero.  If the


          eigenvalue is complex, then WI is SCALE1 times the imaginary


          part of the eigenvalues.  WI will always be non-negative.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAG2S converts a DOUBLE PRECISION matrix, A, to a SINGLE


 PRECISION matrix, SA.


 RMAX is the overflow for the SINGLE PRECISION arithmetic


 DLAG2S checks that all the entries of A are between -RMAX and


 RMAX. If not the conversion is aborted and a flag is raised.


 This is an auxiliary routine so there is no argument checking.

M

          M is INTEGER


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

N

          N is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


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

LDA

          LDA is INTEGER


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

SA

          SA is REAL array, dimension (LDSA,N)


          On exit, if INFO=0, the M-by-N coefficient matrix SA; if


          INFO>0, the content of SA is unspecified.

LDSA

          LDSA is INTEGER


          The leading dimension of the array SA.  LDSA >= max(1,M).

INFO

          INFO is INTEGER


          = 0:  successful exit.


          = 1:  an entry of the matrix A is greater than the SINGLE


                PRECISION overflow threshold, in this case, the content


                of SA in exit is unspecified.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such


 that if ( UPPER ) then


           U**T *A*Q = U**T *( A1 A2 )*Q = ( x  0  )


                             ( 0  A3 )     ( x  x  )


 and


           V**T*B*Q = V**T *( B1 B2 )*Q = ( x  0  )


                            ( 0  B3 )     ( x  x  )


 or if ( .NOT.UPPER ) then


           U**T *A*Q = U**T *( A1 0  )*Q = ( x  x  )


                             ( A2 A3 )     ( 0  x  )


 and


           V**T*B*Q = V**T*( B1 0  )*Q = ( x  x  )


                           ( B2 B3 )     ( 0  x  )


 The rows of the transformed A and B are parallel, where


   U = (  CSU  SNU ), V = (  CSV SNV ), Q = (  CSQ   SNQ )


       ( -SNU  CSU )      ( -SNV CSV )      ( -SNQ   CSQ )


 Z**T denotes the transpose of Z.

UPPER

          UPPER is LOGICAL


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


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

A1

          A1 is DOUBLE PRECISION

A2

          A2 is DOUBLE PRECISION

A3

          A3 is DOUBLE PRECISION


          On entry, A1, A2 and A3 are elements of the input 2-by-2


          upper (lower) triangular matrix A.

B1

          B1 is DOUBLE PRECISION

B2

          B2 is DOUBLE PRECISION

B3

          B3 is DOUBLE PRECISION


          On entry, B1, B2 and B3 are elements of the input 2-by-2


          upper (lower) triangular matrix B.

CSU

          CSU is DOUBLE PRECISION

SNU

          SNU is DOUBLE PRECISION


          The desired orthogonal matrix U.

CSV

          CSV is DOUBLE PRECISION

SNV

          SNV is DOUBLE PRECISION


          The desired orthogonal matrix V.

CSQ

          CSQ is DOUBLE PRECISION

SNQ

          SNQ is DOUBLE PRECISION


          The desired orthogonal matrix Q.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAGTM performs a matrix-vector product of the form


    B := alpha * A * X + beta * B


 where A is a tridiagonal matrix of order N, B and X are N by NRHS


 matrices, and alpha and beta are real scalars, each of which may be


 0., 1., or -1.

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


          = 'C':  Conjugate transpose = Transpose

N

          N is INTEGER


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

NRHS

          NRHS is INTEGER


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


          of the matrices X and B.

ALPHA

          ALPHA is DOUBLE PRECISION


          The scalar alpha.  ALPHA must be 0., 1., or -1.; otherwise,


          it is assumed to be 0.

DL

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


          The (n-1) sub-diagonal elements of T.

D

          D is DOUBLE PRECISION array, dimension (N)


          The diagonal elements of T.

DU

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


          The (n-1) super-diagonal elements of T.

X

          X is DOUBLE PRECISION array, dimension (LDX,NRHS)


          The N by NRHS matrix X.

LDX

          LDX is INTEGER


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

BETA

          BETA is DOUBLE PRECISION


          The scalar beta.  BETA must be 0., 1., or -1.; otherwise,


          it is assumed to be 1.

B

          B is DOUBLE PRECISION array, dimension (LDB,NRHS)


          On entry, the N by NRHS matrix B.


          On exit, B is overwritten by the matrix expression


          B := alpha * A * X + beta * B.

LDB

          LDB is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAGV2 computes the Generalized Schur factorization of a real 2-by-2


 matrix pencil (A,B) where B is upper triangular. This routine


 computes orthogonal (rotation) matrices given by CSL, SNL and CSR,


 SNR such that


 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0


    types), then


    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]


    [  0  a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]


    [ b11 b12 ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]


    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ],


 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues,


    then


    [ a11 a12 ] := [  CSL  SNL ] [ a11 a12 ] [  CSR -SNR ]


    [ a21 a22 ]    [ -SNL  CSL ] [ a21 a22 ] [  SNR  CSR ]


    [ b11  0  ] := [  CSL  SNL ] [ b11 b12 ] [  CSR -SNR ]


    [  0  b22 ]    [ -SNL  CSL ] [  0  b22 ] [  SNR  CSR ]


    where b11 >= b22 > 0.

A

          A is DOUBLE PRECISION array, dimension (LDA, 2)


          On entry, the 2 x 2 matrix A.


          On exit, A is overwritten by the ``A-part'' of the


          generalized Schur form.

LDA

          LDA is INTEGER


          THe leading dimension of the array A.  LDA >= 2.

B

          B is DOUBLE PRECISION array, dimension (LDB, 2)


          On entry, the upper triangular 2 x 2 matrix B.


          On exit, B is overwritten by the ``B-part'' of the


          generalized Schur form.

LDB

          LDB is INTEGER


          THe leading dimension of the array B.  LDB >= 2.

ALPHAR

          ALPHAR is DOUBLE PRECISION array, dimension (2)

ALPHAI

          ALPHAI is DOUBLE PRECISION array, dimension (2)

BETA

          BETA is DOUBLE PRECISION array, dimension (2)


          (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the


          pencil (A,B), k=1,2, i = sqrt(-1).  Note that BETA(k) may


          be zero.

CSL

          CSL is DOUBLE PRECISION


          The cosine of the left rotation matrix.

SNL

          SNL is DOUBLE PRECISION


          The sine of the left rotation matrix.

CSR

          CSR is DOUBLE PRECISION


          The cosine of the right rotation matrix.

SNR

          SNR is DOUBLE PRECISION


          The sine of the right rotation matrix.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

    DLAHQR is an auxiliary routine called by DHSEQR to update the


    eigenvalues and Schur decomposition already computed by DHSEQR, by


    dealing with the Hessenberg submatrix in rows and columns ILO to


    IHI.

WANTT

          WANTT is LOGICAL


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


          = .FALSE.: only eigenvalues are required.

WANTZ

          WANTZ is LOGICAL


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


          = .FALSE.: Schur vectors are not required.

N

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


          It is assumed that H is already upper quasi-triangular in


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


          ILO = 1). DLAHQR works primarily with the Hessenberg


          submatrix in rows and columns ILO to IHI, but applies


          transformations to all of H if WANTT is .TRUE..


          1 <= ILO <= max(1,IHI); IHI <= N.

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


          On entry, the upper Hessenberg matrix H.


          On exit, if INFO is zero and if WANTT is .TRUE., H is upper


          quasi-triangular in rows and columns ILO:IHI, with any


          2-by-2 diagonal blocks in standard form. If INFO is zero


          and WANTT is .FALSE., the contents of H are unspecified on


          exit.  The output state of H if INFO is nonzero is given


          below under the description of INFO.

LDH

          LDH is INTEGER


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

WR

          WR is DOUBLE PRECISION array, dimension (N)

WI

          WI is DOUBLE PRECISION array, dimension (N)


          The real and imaginary parts, respectively, of the computed


          eigenvalues ILO to IHI are stored in the corresponding


          elements of WR and WI. If two eigenvalues are computed as a


          complex conjugate pair, they are stored in consecutive


          elements of WR and WI, say the i-th and (i+1)th, with


          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the


          eigenvalues are stored in the same order as on the diagonal


          of the Schur form returned in H, with WR(i) = H(i,i), and, if


          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,


          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).

ILOZ

          ILOZ is INTEGER

IHIZ

          IHIZ is INTEGER


          Specify the rows of Z to which transformations must be


          applied if WANTZ is .TRUE..


          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ,N)


          If WANTZ is .TRUE., on entry Z must contain the current


          matrix Z of transformations accumulated by DHSEQR, and on


          exit Z has been updated; transformations are applied only to


          the submatrix Z(ILOZ:IHIZ,ILO:IHI).


          If WANTZ is .FALSE., Z is not referenced.

LDZ

          LDZ is INTEGER


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

INFO

          INFO is INTEGER


           = 0:  successful exit


           > 0:  If INFO = i, DLAHQR failed to compute all the


                  eigenvalues ILO to IHI in a total of 30 iterations


                  per eigenvalue; elements i+1:ihi of WR and WI


                  contain those eigenvalues which have been


                  successfully computed.


                  If INFO > 0 and WANTT is .FALSE., then on exit,


                  the remaining unconverged eigenvalues are the


                  eigenvalues of the upper Hessenberg matrix rows


                  and columns ILO through INFO of the final, output


                  value of H.


                  If INFO > 0 and WANTT is .TRUE., then on exit


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


                  where U is an orthogonal matrix.    The final


                  value of H is upper Hessenberg and triangular in


                  rows and columns INFO+1 through IHI.


                  If INFO > 0 and WANTZ is .TRUE., then on exit


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


                  where U is the orthogonal matrix in (*)


                  (regardless of the value of WANTT.)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

     02-96 Based on modifications by


     David Day, Sandia National Laboratory, USA


     12-04 Further modifications by


     Ralph Byers, University of Kansas, USA


     This is a modified version of DLAHQR from LAPACK version 3.0.


     It is (1) more robust against overflow and underflow and


     (2) adopts the more conservative Ahues & Tisseur stopping


     criterion (LAWN 122, 1997).

 DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)


 matrix A so that elements below the k-th subdiagonal are zero. The


 reduction is performed by an orthogonal similarity transformation


 Q**T * A * Q. The routine returns the matrices V and T which determine


 Q as a block reflector I - V*T*V**T, and also the matrix Y = A * V * T.


 This is an auxiliary routine called by DGEHRD.

N

          N is INTEGER


          The order of the matrix A.

K

          K is INTEGER


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


          subdiagonal in the first NB columns are reduced to zero.


          K < N.

NB

          NB is INTEGER


          The number of columns to be reduced.

A

          A is DOUBLE PRECISION array, dimension (LDA,N-K+1)


          On entry, the n-by-(n-k+1) general matrix A.


          On exit, the elements on and above the k-th subdiagonal in


          the first NB columns are overwritten with the corresponding


          elements of the reduced matrix; the elements below the k-th


          subdiagonal, with the array TAU, represent the matrix Q as a


          product of elementary reflectors. The other columns of A are


          unchanged. See Further Details.

LDA

          LDA is INTEGER


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

TAU

          TAU is DOUBLE PRECISION array, dimension (NB)


          The scalar factors of the elementary reflectors. See Further


          Details.

T

          T is DOUBLE PRECISION array, dimension (LDT,NB)


          The upper triangular matrix T.

LDT

          LDT is INTEGER


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

Y

          Y is DOUBLE PRECISION array, dimension (LDY,NB)


          The n-by-nb matrix Y.

LDY

          LDY is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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


  Each H(i) has the form


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


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


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


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


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


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


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


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


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


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


     ( a   a   a   a   a )


     ( a   a   a   a   a )


     ( a   a   a   a   a )


     ( h   h   a   a   a )


     ( v1  h   a   a   a )


     ( v1  v2  a   a   a )


     ( v1  v2  a   a   a )


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


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


  element of the vector defining H(i).


  This subroutine is a slight modification of LAPACK-3.0's DLAHRD


  incorporating improvements proposed by Quintana-Orti and Van de


  Gejin. Note that the entries of A(1:K,2:NB) differ from those


  returned by the original LAPACK-3.0's DLAHRD routine. (This


  subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)

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

 DLAIC1 applies one step of incremental condition estimation in


 its simplest version:


 Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j


 lower triangular matrix L, such that


          twonorm(L*x) = sest


 Then DLAIC1 computes sestpr, s, c such that


 the vector


                 [ s*x ]


          xhat = [  c  ]


 is an approximate singular vector of


                 [ L       0  ]


          Lhat = [ w**T gamma ]


 in the sense that


          twonorm(Lhat*xhat) = sestpr.


 Depending on JOB, an estimate for the largest or smallest singular


 value is computed.


 Note that [s c]**T and sestpr**2 is an eigenpair of the system


     diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]


                                           [ gamma ]


 where  alpha =  x**T*w.

JOB

          JOB is INTEGER


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


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

J

          J is INTEGER


          Length of X and W

X

          X is DOUBLE PRECISION array, dimension (J)


          The j-vector x.

SEST

          SEST is DOUBLE PRECISION


          Estimated singular value of j by j matrix L

W

          W is DOUBLE PRECISION array, dimension (J)


          The j-vector w.

GAMMA

          GAMMA is DOUBLE PRECISION


          The diagonal element gamma.

SESTPR

          SESTPR is DOUBLE PRECISION


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

S

          S is DOUBLE PRECISION


          Sine needed in forming xhat.

C

          C is DOUBLE PRECISION


          Cosine needed in forming xhat.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLALN2 solves a system of the form  (ca A - w D ) X = s B


 or (ca A**T - w D) X = s B   with possible scaling ('s') and


 perturbation of A.  (A**T means A-transpose.)


 A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA


 real diagonal matrix, w is a real or complex value, and X and B are


 NA x 1 matrices -- real if w is real, complex if w is complex.  NA


 may be 1 or 2.


 If w is complex, X and B are represented as NA x 2 matrices,


 the first column of each being the real part and the second


 being the imaginary part.


 's' is a scaling factor (<= 1), computed by DLALN2, which is


 so chosen that X can be computed without overflow.  X is further


 scaled if necessary to assure that norm(ca A - w D)*norm(X) is less


 than overflow.


 If both singular values of (ca A - w D) are less than SMIN,


 SMIN*identity will be used instead of (ca A - w D).  If only one


 singular value is less than SMIN, one element of (ca A - w D) will be


 perturbed enough to make the smallest singular value roughly SMIN.


 If both singular values are at least SMIN, (ca A - w D) will not be


 perturbed.  In any case, the perturbation will be at most some small


 multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values


 are computed by infinity-norm approximations, and thus will only be


 correct to a factor of 2 or so.


 Note: all input quantities are assumed to be smaller than overflow


 by a reasonable factor.  (See BIGNUM.)

LTRANS

          LTRANS is LOGICAL


          =.TRUE.:  A-transpose will be used.


          =.FALSE.: A will be used (not transposed.)

NA

          NA is INTEGER


          The size of the matrix A.  It may (only) be 1 or 2.

NW

          NW is INTEGER


          1 if 'w' is real, 2 if 'w' is complex.  It may only be 1


          or 2.

SMIN

          SMIN is DOUBLE PRECISION


          The desired lower bound on the singular values of A.  This


          should be a safe distance away from underflow or overflow,


          say, between (underflow/machine precision) and  (machine


          precision * overflow ).  (See BIGNUM and ULP.)

CA

          CA is DOUBLE PRECISION


          The coefficient c, which A is multiplied by.

A

          A is DOUBLE PRECISION array, dimension (LDA,NA)


          The NA x NA matrix A.

LDA

          LDA is INTEGER


          The leading dimension of A.  It must be at least NA.

D1

          D1 is DOUBLE PRECISION


          The 1,1 element in the diagonal matrix D.

D2

          D2 is DOUBLE PRECISION


          The 2,2 element in the diagonal matrix D.  Not used if NA=1.

B

          B is DOUBLE PRECISION array, dimension (LDB,NW)


          The NA x NW matrix B (right-hand side).  If NW=2 ('w' is


          complex), column 1 contains the real part of B and column 2


          contains the imaginary part.

LDB

          LDB is INTEGER


          The leading dimension of B.  It must be at least NA.

WR

          WR is DOUBLE PRECISION


          The real part of the scalar 'w'.

WI

          WI is DOUBLE PRECISION


          The imaginary part of the scalar 'w'.  Not used if NW=1.

X

          X is DOUBLE PRECISION array, dimension (LDX,NW)


          The NA x NW matrix X (unknowns), as computed by DLALN2.


          If NW=2 ('w' is complex), on exit, column 1 will contain


          the real part of X and column 2 will contain the imaginary


          part.

LDX

          LDX is INTEGER


          The leading dimension of X.  It must be at least NA.

SCALE

          SCALE is DOUBLE PRECISION


          The scale factor that B must be multiplied by to insure


          that overflow does not occur when computing X.  Thus,


          (ca A - w D) X  will be SCALE*B, not B (ignoring


          perturbations of A.)  It will be at most 1.

XNORM

          XNORM is DOUBLE PRECISION


          The infinity-norm of X, when X is regarded as an NA x NW


          real matrix.

INFO

          INFO is INTEGER


          An error flag.  It will be set to zero if no error occurs,


          a negative number if an argument is in error, or a positive


          number if  ca A - w D  had to be perturbed.


          The possible values are:


          = 0: No error occurred, and (ca A - w D) did not have to be


                 perturbed.


          = 1: (ca A - w D) had to be perturbed to make its smallest


               (or only) singular value greater than SMIN.


          NOTE: In the interests of speed, this routine does not


                check the inputs for errors.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANGT  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the  element of  largest absolute value  of a


 real tridiagonal matrix A.

DLANGT

    DLANGT = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANGT as described


          above.

N

          N is INTEGER


          The order of the matrix A.  N >= 0.  When N = 0, DLANGT is


          set to zero.

DL

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


          The (n-1) sub-diagonal elements of A.

D

          D is DOUBLE PRECISION array, dimension (N)


          The diagonal elements of A.

DU

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


          The (n-1) super-diagonal elements of A.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANHS  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the  element of  largest absolute value  of a


 Hessenberg matrix A.

DLANHS

    DLANHS = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANHS as described


          above.

N

          N is INTEGER


          The order of the matrix A.  N >= 0.  When N = 0, DLANHS is


          set to zero.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          The n by n upper Hessenberg matrix A; the part of A below the


          first sub-diagonal is not referenced.

LDA

          LDA is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


          where LWORK >= N when NORM = 'I'; otherwise, WORK is not


          referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANSB  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the element of  largest absolute value  of an


 n by n symmetric band matrix A,  with k super-diagonals.

DLANSB

    DLANSB = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANSB as described


          above.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          band matrix A is supplied.


          = 'U':  Upper triangular part is supplied


          = 'L':  Lower triangular part is supplied

N

          N is INTEGER


          The order of the matrix A.  N >= 0.  When N = 0, DLANSB is


          set to zero.

K

          K is INTEGER


          The number of super-diagonals or sub-diagonals of the


          band matrix A.  K >= 0.

AB

          AB is DOUBLE PRECISION array, dimension (LDAB,N)


          The upper or lower triangle of the symmetric band matrix A,


          stored in the first K+1 rows of AB.  The j-th column of A is


          stored in the j-th column of the array AB as follows:


          if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;


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

LDAB

          LDAB is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,


          WORK is not referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANSP  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the  element of  largest absolute value  of a


 real symmetric matrix A,  supplied in packed form.

DLANSP

    DLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANSP as described


          above.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          symmetric matrix A is supplied.


          = 'U':  Upper triangular part of A is supplied


          = 'L':  Lower triangular part of A is supplied

N

          N is INTEGER


          The order of the matrix A.  N >= 0.  When N = 0, DLANSP is


          set to zero.

AP

          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)


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


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


          in the array AP as follows:


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


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

WORK

          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,


          WORK is not referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANTB  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the element of  largest absolute value  of an


 n by n triangular band matrix A,  with ( k + 1 ) diagonals.

DLANTB

    DLANTB = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANTB as described


          above.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

N

          N is INTEGER


          The order of the matrix A.  N >= 0.  When N = 0, DLANTB is


          set to zero.

K

          K is INTEGER


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


          or the number of sub-diagonals of the matrix A if UPLO = 'L'.


          K >= 0.

AB

          AB is DOUBLE PRECISION array, dimension (LDAB,N)


          The upper or lower triangular band matrix A, stored in the


          first k+1 rows of AB.  The j-th column of A is stored


          in the j-th column of the array AB as follows:


          if UPLO = 'U', AB(k+1+i-j,j) = A(i,j) for max(1,j-k)<=i<=j;


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


          Note that when DIAG = 'U', the elements of the array AB


          corresponding to the diagonal elements of the matrix A are


          not referenced, but are assumed to be one.

LDAB

          LDAB is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


          where LWORK >= N when NORM = 'I'; otherwise, WORK is not


          referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANTP  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the  element of  largest absolute value  of a


 triangular matrix A, supplied in packed form.

DLANTP

    DLANTP = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANTP as described


          above.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

N

          N is INTEGER


          The order of the matrix A.  N >= 0.  When N = 0, DLANTP is


          set to zero.

AP

          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)


          The upper or lower triangular matrix A, packed columnwise in


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


          AP as follows:


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


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


          Note that when DIAG = 'U', the elements of the array AP


          corresponding to the diagonal elements of the matrix A are


          not referenced, but are assumed to be one.

WORK

          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


          where LWORK >= N when NORM = 'I'; otherwise, WORK is not


          referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANTR  returns the value of the one norm,  or the Frobenius norm, or


 the  infinity norm,  or the  element of  largest absolute value  of a


 trapezoidal or triangular matrix A.

DLANTR

    DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm'


             (


             ( norm1(A),         NORM = '1', 'O' or 'o'


             (


             ( normI(A),         NORM = 'I' or 'i'


             (


             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'


 where  norm1  denotes the  one norm of a matrix (maximum column sum),


 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and


 normF  denotes the  Frobenius norm of a matrix (square root of sum of


 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.

NORM

          NORM is CHARACTER*1


          Specifies the value to be returned in DLANTR as described


          above.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower trapezoidal.


          = 'U':  Upper trapezoidal


          = 'L':  Lower trapezoidal


          Note that A is triangular instead of trapezoidal if M = N.

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A has unit diagonal.


          = 'N':  Non-unit diagonal


          = 'U':  Unit diagonal

M

          M is INTEGER


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


          UPLO = 'U', M <= N.  When M = 0, DLANTR is set to zero.

N

          N is INTEGER


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


          UPLO = 'L', N <= M.  When N = 0, DLANTR is set to zero.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          The trapezoidal matrix A (A is triangular if M = N).


          If UPLO = 'U', the leading m by n upper trapezoidal part of


          the array A contains the upper trapezoidal matrix, and the


          strictly lower triangular part of A is not referenced.


          If UPLO = 'L', the leading m by n lower trapezoidal part of


          the array A contains the lower trapezoidal matrix, and the


          strictly upper triangular part of A is not referenced.  Note


          that when DIAG = 'U', the diagonal elements of A are not


          referenced and are assumed to be one.

LDA

          LDA is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


          where LWORK >= M when NORM = 'I'; otherwise, WORK is not


          referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric


 matrix in standard form:


      [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]


      [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]


 where either


 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or


 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex


 conjugate eigenvalues.

A

          A is DOUBLE PRECISION

B

          B is DOUBLE PRECISION

C

          C is DOUBLE PRECISION

D

          D is DOUBLE PRECISION


          On entry, the elements of the input matrix.


          On exit, they are overwritten by the elements of the


          standardised Schur form.

RT1R

          RT1R is DOUBLE PRECISION

RT1I

          RT1I is DOUBLE PRECISION

RT2R

          RT2R is DOUBLE PRECISION

RT2I

          RT2I is DOUBLE PRECISION


          The real and imaginary parts of the eigenvalues. If the


          eigenvalues are a complex conjugate pair, RT1I > 0.

CS

          CS is DOUBLE PRECISION

SN

          SN is DOUBLE PRECISION


          Parameters of the rotation matrix.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Modified by V. Sima, Research Institute for Informatics, Bucharest,


  Romania, to reduce the risk of cancellation errors,


  when computing real eigenvalues, and to ensure, if possible, that


  abs(RT1R) >= abs(RT2R).

 Given two column vectors X and Y, let


                      A = ( X Y ).


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


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


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


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

N

          N is INTEGER


          The length of the vectors X and Y.

X

          X is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCX)


          On entry, X contains the N-vector X.


          On exit, X is overwritten.

INCX

          INCX is INTEGER


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

Y

          Y is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCY)


          On entry, Y contains the N-vector Y.


          On exit, Y is overwritten.

INCY

          INCY is INTEGER


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

SSMIN

          SSMIN is DOUBLE PRECISION


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAPMR rearranges the rows of the M by N matrix X as specified


 by the permutation K(1),K(2),...,K(M) of the integers 1,...,M.


 If FORWRD = .TRUE.,  forward permutation:


      X(K(I),*) is moved X(I,*) for I = 1,2,...,M.


 If FORWRD = .FALSE., backward permutation:


      X(I,*) is moved to X(K(I),*) for I = 1,2,...,M.

FORWRD

          FORWRD is LOGICAL


          = .TRUE., forward permutation


          = .FALSE., backward permutation

M

          M is INTEGER


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

N

          N is INTEGER


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

X

          X is DOUBLE PRECISION array, dimension (LDX,N)


          On entry, the M by N matrix X.


          On exit, X contains the permuted matrix X.

LDX

          LDX is INTEGER


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

K

          K is INTEGER array, dimension (M)


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


          internal workspace, but reset to its original value on


          output.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAPMT rearranges the columns of the M by N matrix X as specified


 by the permutation K(1),K(2),...,K(N) of the integers 1,...,N.


 If FORWRD = .TRUE.,  forward permutation:


      X(*,K(J)) is moved X(*,J) for J = 1,2,...,N.


 If FORWRD = .FALSE., backward permutation:


      X(*,J) is moved to X(*,K(J)) for J = 1,2,...,N.

FORWRD

          FORWRD is LOGICAL


          = .TRUE., forward permutation


          = .FALSE., backward permutation

M

          M is INTEGER


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

N

          N is INTEGER


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

X

          X is DOUBLE PRECISION array, dimension (LDX,N)


          On entry, the M by N matrix X.


          On exit, X contains the permuted matrix X.

LDX

          LDX is INTEGER


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

K

          K is INTEGER array, dimension (N)


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


          internal workspace, but reset to its original value on


          output.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAQP2 computes a QR factorization with column pivoting of


 the block A(OFFSET+1:M,1:N).


 The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.

M

          M is INTEGER


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

N

          N is INTEGER


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

OFFSET

          OFFSET is INTEGER


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


          but no factorized. OFFSET >= 0.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


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


          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is


          the triangular factor obtained; the elements in block


          A(OFFSET+1:M,1:N) below the diagonal, together with the


          array TAU, represent the orthogonal matrix Q as a product of


          elementary reflectors. Block A(1:OFFSET,1:N) has been


          accordingly pivoted, but no factorized.

LDA

          LDA is INTEGER


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

JPVT

          JPVT is INTEGER array, dimension (N)


          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted


          to the front of A*P (a leading column); if JPVT(i) = 0,


          the i-th column of A is a free column.


          On exit, if JPVT(i) = k, then the i-th column of A*P


          was the k-th column of A.

TAU

          TAU is DOUBLE PRECISION array, dimension (min(M,N))


          The scalar factors of the elementary reflectors.

VN1

          VN1 is DOUBLE PRECISION array, dimension (N)


          The vector with the partial column norms.

VN2

          VN2 is DOUBLE PRECISION array, dimension (N)


          The vector with the exact column norms.

WORK

          WORK is DOUBLE PRECISION array, dimension (N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

LAPACK Working Note 176

 DLAQPS computes a step of QR factorization with column pivoting


 of a real M-by-N matrix A by using Blas-3.  It tries to factorize


 NB columns from A starting from the row OFFSET+1, and updates all


 of the matrix with Blas-3 xGEMM.


 In some cases, due to catastrophic cancellations, it cannot


 factorize NB columns.  Hence, the actual number of factorized


 columns is returned in KB.


 Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.

M

          M is INTEGER


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

N

          N is INTEGER


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

OFFSET

          OFFSET is INTEGER


          The number of rows of A that have been factorized in


          previous steps.

NB

          NB is INTEGER


          The number of columns to factorize.

KB

          KB is INTEGER


          The number of columns actually factorized.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


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


          On exit, block A(OFFSET+1:M,1:KB) is the triangular


          factor obtained and block A(1:OFFSET,1:N) has been


          accordingly pivoted, but no factorized.


          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has


          been updated.

LDA

          LDA is INTEGER


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

JPVT

          JPVT is INTEGER array, dimension (N)


          JPVT(I) = K <==> Column K of the full matrix A has been


          permuted into position I in AP.

TAU

          TAU is DOUBLE PRECISION array, dimension (KB)


          The scalar factors of the elementary reflectors.

VN1

          VN1 is DOUBLE PRECISION array, dimension (N)


          The vector with the partial column norms.

VN2

          VN2 is DOUBLE PRECISION array, dimension (N)


          The vector with the exact column norms.

AUXV

          AUXV is DOUBLE PRECISION array, dimension (NB)


          Auxiliary vector.

F

          F is DOUBLE PRECISION array, dimension (LDF,NB)


          Matrix F**T = L*Y**T*A.

LDF

          LDF is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

LAPACK Working Note 176

    DLAQR0 computes the eigenvalues of a Hessenberg matrix H


    and, optionally, the matrices T and Z from the Schur decomposition


    H = Z T Z**T, where T is an upper quasi-triangular matrix (the


    Schur form), and Z is the orthogonal matrix of Schur vectors.


    Optionally Z may be postmultiplied into an input orthogonal


    matrix Q so that this routine can give the Schur factorization


    of a matrix A which has been reduced to the Hessenberg form H


    by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.

WANTT

          WANTT is LOGICAL


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


          = .FALSE.: only eigenvalues are required.

WANTZ

          WANTZ is LOGICAL


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


          = .FALSE.: Schur vectors are not required.

N

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


           It is assumed that H is already upper triangular in rows


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


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


           previous call to DGEBAL, and then passed to DGEHRD when the


           matrix output by DGEBAL is reduced to Hessenberg form.


           Otherwise, ILO and IHI should be set to 1 and N,


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


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

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


           On entry, the upper Hessenberg matrix H.


           On exit, if INFO = 0 and WANTT is .TRUE., then H contains


           the upper quasi-triangular matrix T from the Schur


           decomposition (the Schur form); 2-by-2 diagonal blocks


           (corresponding to complex conjugate pairs of eigenvalues)


           are returned in standard form, with H(i,i) = H(i+1,i+1)


           and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is


           .FALSE., then the contents of H are unspecified on exit.


           (The output value of H when INFO > 0 is given under the


           description of INFO below.)


           This subroutine may explicitly set H(i,j) = 0 for i > j and


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

LDH

          LDH is INTEGER


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

WR

          WR is DOUBLE PRECISION array, dimension (IHI)

WI

          WI is DOUBLE PRECISION array, dimension (IHI)


           The real and imaginary parts, respectively, of the computed


           eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)


           and WI(ILO:IHI). If two eigenvalues are computed as a


           complex conjugate pair, they are stored in consecutive


           elements of WR and WI, say the i-th and (i+1)th, with


           WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then


           the eigenvalues are stored in the same order as on the


           diagonal of the Schur form returned in H, with


           WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal


           block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and


           WI(i+1) = -WI(i).

ILOZ

          ILOZ is INTEGER

IHIZ

          IHIZ is INTEGER


           Specify the rows of Z to which transformations must be


           applied if WANTZ is .TRUE..


           1 <= ILOZ <= ILO; IHI <= IHIZ <= N.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ,IHI)


           If WANTZ is .FALSE., then Z is not referenced.


           If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is


           replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the


           orthogonal Schur factor of H(ILO:IHI,ILO:IHI).


           (The output value of Z when INFO > 0 is given under


           the description of INFO below.)

LDZ

          LDZ is INTEGER


           The leading dimension of the array Z.  if WANTZ is .TRUE.


           then LDZ >= MAX(1,IHIZ).  Otherwise, LDZ >= 1.

WORK

          WORK is DOUBLE PRECISION array, dimension LWORK


           On exit, if LWORK = -1, WORK(1) returns an estimate of


           the optimal value for LWORK.

LWORK

          LWORK is INTEGER


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


           is sufficient, but LWORK typically as large as 6*N may


           be required for optimal performance.  A workspace query


           to determine the optimal workspace size is recommended.


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


           In this case, DLAQR0 checks the input parameters and


           estimates the optimal workspace size for the given


           values of N, ILO and IHI.  The estimate is returned


           in WORK(1).  No error message related to LWORK is


           issued by XERBLA.  Neither H nor Z are accessed.

INFO

          INFO is INTEGER


             = 0:  successful exit


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


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


                and WI contain those eigenvalues which have been


                successfully computed.  (Failures are rare.)


                If INFO > 0 and WANT is .FALSE., then on exit,


                the remaining unconverged eigenvalues are the eigen-


                values of the upper Hessenberg matrix rows and


                columns ILO through INFO of the final, output


                value of H.


                If INFO > 0 and WANTT is .TRUE., then on exit


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


                where U is an orthogonal matrix.  The final


                value of H is upper Hessenberg and quasi-triangular


                in rows and columns INFO+1 through IHI.


                If INFO > 0 and WANTZ is .TRUE., then on exit


                  (final value of Z(ILO:IHI,ILOZ:IHIZ)


                   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U


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


                less of the value of WANTT.)


                If INFO > 0 and WANTZ is .FALSE., then Z is not


                accessed.

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

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


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


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


  929--947, 2002.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

      Given a 2-by-2 or 3-by-3 matrix H, DLAQR1 sets v to a


      scalar multiple of the first column of the product


      (*)  K = (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)


      scaling to avoid overflows and most underflows. It


      is assumed that either


              1) sr1 = sr2 and si1 = -si2


          or


              2) si1 = si2 = 0.


      This is useful for starting double implicit shift bulges


      in the QR algorithm.

N

          N is INTEGER


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

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


              The 2-by-2 or 3-by-3 matrix H in (*).

LDH

          LDH is INTEGER


              The leading dimension of H as declared in


              the calling procedure.  LDH >= N

SR1

          SR1 is DOUBLE PRECISION

SI1

          SI1 is DOUBLE PRECISION

SR2

          SR2 is DOUBLE PRECISION

SI2

          SI2 is DOUBLE PRECISION


              The shifts in (*).

V

          V is DOUBLE PRECISION array, dimension (N)


              A scalar multiple of the first column of the


              matrix K in (*).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

    DLAQR2 is identical to DLAQR3 except that it avoids


    recursion by calling DLAHQR instead of DLAQR4.


    Aggressive early deflation:


    This subroutine accepts as input an upper Hessenberg matrix


    H and performs an orthogonal similarity transformation


    designed to detect and deflate fully converged eigenvalues from


    a trailing principal submatrix.  On output H has been over-


    written by a new Hessenberg matrix that is a perturbation of


    an orthogonal similarity transformation of H.  It is to be


    hoped that the final version of H has many zero subdiagonal


    entries.

WANTT

          WANTT is LOGICAL


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


          so that the quasi-triangular Schur factor may be


          computed (in cooperation with the calling subroutine).


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


          the eigenvalues.

WANTZ

          WANTZ is LOGICAL


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


          so that the orthogonal Schur factor may be computed


          (in cooperation with the calling subroutine).


          If .FALSE., then Z is not referenced.

N

          N is INTEGER


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


          order of the orthogonal matrix Z.

KTOP

          KTOP is INTEGER


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


          KBOT and KTOP together determine an isolated block


          along the diagonal of the Hessenberg matrix.

KBOT

          KBOT is INTEGER


          It is assumed without a check that either


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


          determine an isolated block along the diagonal of the


          Hessenberg matrix.

NW

          NW is INTEGER


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

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


          On input the initial N-by-N section of H stores the


          Hessenberg matrix undergoing aggressive early deflation.


          On output H has been transformed by an orthogonal


          similarity transformation, perturbed, and the returned


          to Hessenberg form that (it is to be hoped) has some


          zero subdiagonal entries.

LDH

          LDH is INTEGER


          Leading dimension of H just as declared in the calling


          subroutine.  N <= LDH

ILOZ

          ILOZ is INTEGER

IHIZ

          IHIZ is INTEGER


          Specify the rows of Z to which transformations must be


          applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ,N)


          IF WANTZ is .TRUE., then on output, the orthogonal


          similarity transformation mentioned above has been


          accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.


          If WANTZ is .FALSE., then Z is unreferenced.

LDZ

          LDZ is INTEGER


          The leading dimension of Z just as declared in the


          calling subroutine.  1 <= LDZ.

NS

          NS is INTEGER


          The number of unconverged (ie approximate) eigenvalues


          returned in SR and SI that may be used as shifts by the


          calling subroutine.

ND

          ND is INTEGER


          The number of converged eigenvalues uncovered by this


          subroutine.

SR

          SR is DOUBLE PRECISION array, dimension (KBOT)

SI

          SI is DOUBLE PRECISION array, dimension (KBOT)


          On output, the real and imaginary parts of approximate


          eigenvalues that may be used for shifts are stored in


          SR(KBOT-ND-NS+1) through SR(KBOT-ND) and


          SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.


          The real and imaginary parts of converged eigenvalues


          are stored in SR(KBOT-ND+1) through SR(KBOT) and


          SI(KBOT-ND+1) through SI(KBOT), respectively.

V

          V is DOUBLE PRECISION array, dimension (LDV,NW)


          An NW-by-NW work array.

LDV

          LDV is INTEGER


          The leading dimension of V just as declared in the


          calling subroutine.  NW <= LDV

NH

          NH is INTEGER


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

T

          T is DOUBLE PRECISION array, dimension (LDT,NW)

LDT

          LDT is INTEGER


          The leading dimension of T just as declared in the


          calling subroutine.  NW <= LDT

NV

          NV is INTEGER


          The number of rows of work array WV available for


          workspace.  NV >= NW.

WV

          WV is DOUBLE PRECISION array, dimension (LDWV,NW)

LDWV

          LDWV is INTEGER


          The leading dimension of W just as declared in the


          calling subroutine.  NW <= LDV

WORK

          WORK is DOUBLE PRECISION array, dimension (LWORK)


          On exit, WORK(1) is set to an estimate of the optimal value


          of LWORK for the given values of N, NW, KTOP and KBOT.

LWORK

          LWORK is INTEGER


          The dimension of the work array WORK.  LWORK = 2*NW


          suffices, but greater efficiency may result from larger


          values of LWORK.


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


          only estimates the optimal workspace size for the given


          values of N, NW, KTOP and KBOT.  The estimate is returned


          in WORK(1).  No error message related to LWORK is issued


          by XERBLA.  Neither H nor Z are accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

    Aggressive early deflation:


    DLAQR3 accepts as input an upper Hessenberg matrix


    H and performs an orthogonal similarity transformation


    designed to detect and deflate fully converged eigenvalues from


    a trailing principal submatrix.  On output H has been over-


    written by a new Hessenberg matrix that is a perturbation of


    an orthogonal similarity transformation of H.  It is to be


    hoped that the final version of H has many zero subdiagonal


    entries.

WANTT

          WANTT is LOGICAL


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


          so that the quasi-triangular Schur factor may be


          computed (in cooperation with the calling subroutine).


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


          the eigenvalues.

WANTZ

          WANTZ is LOGICAL


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


          so that the orthogonal Schur factor may be computed


          (in cooperation with the calling subroutine).


          If .FALSE., then Z is not referenced.

N

          N is INTEGER


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


          order of the orthogonal matrix Z.

KTOP

          KTOP is INTEGER


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


          KBOT and KTOP together determine an isolated block


          along the diagonal of the Hessenberg matrix.

KBOT

          KBOT is INTEGER


          It is assumed without a check that either


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


          determine an isolated block along the diagonal of the


          Hessenberg matrix.

NW

          NW is INTEGER


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

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


          On input the initial N-by-N section of H stores the


          Hessenberg matrix undergoing aggressive early deflation.


          On output H has been transformed by an orthogonal


          similarity transformation, perturbed, and the returned


          to Hessenberg form that (it is to be hoped) has some


          zero subdiagonal entries.

LDH

          LDH is INTEGER


          Leading dimension of H just as declared in the calling


          subroutine.  N <= LDH

ILOZ

          ILOZ is INTEGER

IHIZ

          IHIZ is INTEGER


          Specify the rows of Z to which transformations must be


          applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ,N)


          IF WANTZ is .TRUE., then on output, the orthogonal


          similarity transformation mentioned above has been


          accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.


          If WANTZ is .FALSE., then Z is unreferenced.

LDZ

          LDZ is INTEGER


          The leading dimension of Z just as declared in the


          calling subroutine.  1 <= LDZ.

NS

          NS is INTEGER


          The number of unconverged (ie approximate) eigenvalues


          returned in SR and SI that may be used as shifts by the


          calling subroutine.

ND

          ND is INTEGER


          The number of converged eigenvalues uncovered by this


          subroutine.

SR

          SR is DOUBLE PRECISION array, dimension (KBOT)

SI

          SI is DOUBLE PRECISION array, dimension (KBOT)


          On output, the real and imaginary parts of approximate


          eigenvalues that may be used for shifts are stored in


          SR(KBOT-ND-NS+1) through SR(KBOT-ND) and


          SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.


          The real and imaginary parts of converged eigenvalues


          are stored in SR(KBOT-ND+1) through SR(KBOT) and


          SI(KBOT-ND+1) through SI(KBOT), respectively.

V

          V is DOUBLE PRECISION array, dimension (LDV,NW)


          An NW-by-NW work array.

LDV

          LDV is INTEGER


          The leading dimension of V just as declared in the


          calling subroutine.  NW <= LDV

NH

          NH is INTEGER


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

T

          T is DOUBLE PRECISION array, dimension (LDT,NW)

LDT

          LDT is INTEGER


          The leading dimension of T just as declared in the


          calling subroutine.  NW <= LDT

NV

          NV is INTEGER


          The number of rows of work array WV available for


          workspace.  NV >= NW.

WV

          WV is DOUBLE PRECISION array, dimension (LDWV,NW)

LDWV

          LDWV is INTEGER


          The leading dimension of W just as declared in the


          calling subroutine.  NW <= LDV

WORK

          WORK is DOUBLE PRECISION array, dimension (LWORK)


          On exit, WORK(1) is set to an estimate of the optimal value


          of LWORK for the given values of N, NW, KTOP and KBOT.

LWORK

          LWORK is INTEGER


          The dimension of the work array WORK.  LWORK = 2*NW


          suffices, but greater efficiency may result from larger


          values of LWORK.


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


          only estimates the optimal workspace size for the given


          values of N, NW, KTOP and KBOT.  The estimate is returned


          in WORK(1).  No error message related to LWORK is issued


          by XERBLA.  Neither H nor Z are accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

    DLAQR4 implements one level of recursion for DLAQR0.


    It is a complete implementation of the small bulge multi-shift


    QR algorithm.  It may be called by DLAQR0 and, for large enough


    deflation window size, it may be called by DLAQR3.  This


    subroutine is identical to DLAQR0 except that it calls DLAQR2


    instead of DLAQR3.


    DLAQR4 computes the eigenvalues of a Hessenberg matrix H


    and, optionally, the matrices T and Z from the Schur decomposition


    H = Z T Z**T, where T is an upper quasi-triangular matrix (the


    Schur form), and Z is the orthogonal matrix of Schur vectors.


    Optionally Z may be postmultiplied into an input orthogonal


    matrix Q so that this routine can give the Schur factorization


    of a matrix A which has been reduced to the Hessenberg form H


    by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.

WANTT

          WANTT is LOGICAL


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


          = .FALSE.: only eigenvalues are required.

WANTZ

          WANTZ is LOGICAL


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


          = .FALSE.: Schur vectors are not required.

N

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


           It is assumed that H is already upper triangular in rows


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


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


           previous call to DGEBAL, and then passed to DGEHRD when the


           matrix output by DGEBAL is reduced to Hessenberg form.


           Otherwise, ILO and IHI should be set to 1 and N,


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


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

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


           On entry, the upper Hessenberg matrix H.


           On exit, if INFO = 0 and WANTT is .TRUE., then H contains


           the upper quasi-triangular matrix T from the Schur


           decomposition (the Schur form); 2-by-2 diagonal blocks


           (corresponding to complex conjugate pairs of eigenvalues)


           are returned in standard form, with H(i,i) = H(i+1,i+1)


           and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is


           .FALSE., then the contents of H are unspecified on exit.


           (The output value of H when INFO > 0 is given under the


           description of INFO below.)


           This subroutine may explicitly set H(i,j) = 0 for i > j and


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

LDH

          LDH is INTEGER


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

WR

          WR is DOUBLE PRECISION array, dimension (IHI)

WI

          WI is DOUBLE PRECISION array, dimension (IHI)


           The real and imaginary parts, respectively, of the computed


           eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)


           and WI(ILO:IHI). If two eigenvalues are computed as a


           complex conjugate pair, they are stored in consecutive


           elements of WR and WI, say the i-th and (i+1)th, with


           WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then


           the eigenvalues are stored in the same order as on the


           diagonal of the Schur form returned in H, with


           WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal


           block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and


           WI(i+1) = -WI(i).

ILOZ

          ILOZ is INTEGER

IHIZ

          IHIZ is INTEGER


           Specify the rows of Z to which transformations must be


           applied if WANTZ is .TRUE..


           1 <= ILOZ <= ILO; IHI <= IHIZ <= N.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ,IHI)


           If WANTZ is .FALSE., then Z is not referenced.


           If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is


           replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the


           orthogonal Schur factor of H(ILO:IHI,ILO:IHI).


           (The output value of Z when INFO > 0 is given under


           the description of INFO below.)

LDZ

          LDZ is INTEGER


           The leading dimension of the array Z.  if WANTZ is .TRUE.


           then LDZ >= MAX(1,IHIZ).  Otherwise, LDZ >= 1.

WORK

          WORK is DOUBLE PRECISION array, dimension LWORK


           On exit, if LWORK = -1, WORK(1) returns an estimate of


           the optimal value for LWORK.

LWORK

          LWORK is INTEGER


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


           is sufficient, but LWORK typically as large as 6*N may


           be required for optimal performance.  A workspace query


           to determine the optimal workspace size is recommended.


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


           In this case, DLAQR4 checks the input parameters and


           estimates the optimal workspace size for the given


           values of N, ILO and IHI.  The estimate is returned


           in WORK(1).  No error message related to LWORK is


           issued by XERBLA.  Neither H nor Z are accessed.

INFO

          INFO is INTEGER


             = 0:  successful exit


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


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


                and WI contain those eigenvalues which have been


                successfully computed.  (Failures are rare.)


                If INFO > 0 and WANT is .FALSE., then on exit,


                the remaining unconverged eigenvalues are the eigen-


                values of the upper Hessenberg matrix rows and


                columns ILO through INFO of the final, output


                value of H.


                If INFO > 0 and WANTT is .TRUE., then on exit


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


                where U is a orthogonal matrix.  The final


                value of  H is upper Hessenberg and triangular in


                rows and columns INFO+1 through IHI.


                If INFO > 0 and WANTZ is .TRUE., then on exit


                  (final value of Z(ILO:IHI,ILOZ:IHIZ)


                   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U


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


                less of the value of WANTT.)


                If INFO > 0 and WANTZ is .FALSE., then Z is not


                accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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


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


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


  929--947, 2002.

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

    DLAQR5, called by DLAQR0, performs a


    single small-bulge multi-shift QR sweep.

WANTT

          WANTT is LOGICAL


             WANTT = .true. if the quasi-triangular Schur factor


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

WANTZ

          WANTZ is LOGICAL


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


             computed.  WANTZ is set to .false. otherwise.

KACC22

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


             Specifies the computation mode of far-from-diagonal


             orthogonal updates.


        = 0: DLAQR5 does not accumulate reflections and does not


             use matrix-matrix multiply to update far-from-diagonal


             matrix entries.


        = 1: DLAQR5 accumulates reflections and uses matrix-matrix


             multiply to update the far-from-diagonal matrix entries.


        = 2: Same as KACC22 = 1. This option used to enable exploiting


             the 2-by-2 structure during matrix multiplications, but


             this is no longer supported.

N

          N is INTEGER


             N is the order of the Hessenberg matrix H upon which this


             subroutine operates.

KTOP

          KTOP is INTEGER

KBOT

          KBOT is INTEGER


             These are the first and last rows and columns of an


             isolated diagonal block upon which the QR sweep is to be


             applied. It is assumed without a check that


                       either KTOP = 1  or   H(KTOP,KTOP-1) = 0


             and


                       either KBOT = N  or   H(KBOT+1,KBOT) = 0.

NSHFTS

          NSHFTS is INTEGER


             NSHFTS gives the number of simultaneous shifts.  NSHFTS


             must be positive and even.

SR

          SR is DOUBLE PRECISION array, dimension (NSHFTS)

SI

          SI is DOUBLE PRECISION array, dimension (NSHFTS)


             SR contains the real parts and SI contains the imaginary


             parts of the NSHFTS shifts of origin that define the


             multi-shift QR sweep.  On output SR and SI may be


             reordered.

H

          H is DOUBLE PRECISION array, dimension (LDH,N)


             On input H contains a Hessenberg matrix.  On output a


             multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied


             to the isolated diagonal block in rows and columns KTOP


             through KBOT.

LDH

          LDH is INTEGER


             LDH is the leading dimension of H just as declared in the


             calling procedure.  LDH >= MAX(1,N).

ILOZ

          ILOZ is INTEGER

IHIZ

          IHIZ is INTEGER


             Specify the rows of Z to which transformations must be


             applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N

Z

          Z is DOUBLE PRECISION array, dimension (LDZ,IHIZ)


             If WANTZ = .TRUE., then the QR Sweep orthogonal


             similarity transformation is accumulated into


             Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.


             If WANTZ = .FALSE., then Z is unreferenced.

LDZ

          LDZ is INTEGER


             LDA is the leading dimension of Z just as declared in


             the calling procedure. LDZ >= N.

V

          V is DOUBLE PRECISION array, dimension (LDV,NSHFTS/2)

LDV

          LDV is INTEGER


             LDV is the leading dimension of V as declared in the


             calling procedure.  LDV >= 3.

U

          U is DOUBLE PRECISION array, dimension (LDU,2*NSHFTS)

LDU

          LDU is INTEGER


             LDU is the leading dimension of U just as declared in the


             in the calling subroutine.  LDU >= 2*NSHFTS.

NV

          NV is INTEGER


             NV is the number of rows in WV agailable for workspace.


             NV >= 1.

WV

          WV is DOUBLE PRECISION array, dimension (LDWV,2*NSHFTS)

LDWV

          LDWV is INTEGER


             LDWV is the leading dimension of WV as declared in the


             in the calling subroutine.  LDWV >= NV.

NH

          NH is INTEGER


             NH is the number of columns in array WH available for


             workspace. NH >= 1.

WH

          WH is DOUBLE PRECISION array, dimension (LDWH,NH)

LDWH

          LDWH is INTEGER


             Leading dimension of WH just as declared in the


             calling procedure.  LDWH >= 2*NSHFTS.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

 DLAQSB equilibrates a symmetric band matrix A using the scaling


 factors in the vector S.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          symmetric matrix A is stored.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

N

          N is INTEGER


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

KD

          KD is INTEGER


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


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

AB

          AB is DOUBLE PRECISION array, dimension (LDAB,N)


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


          matrix A, stored in the first KD+1 rows of the array.  The


          j-th column of A is stored in the j-th column of the array AB


          as follows:


          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;


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


          On exit, if INFO = 0, the triangular factor U or L from the


          Cholesky factorization A = U**T*U or A = L*L**T of the band


          matrix A, in the same storage format as A.

LDAB

          LDAB is INTEGER


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

S

          S is DOUBLE PRECISION array, dimension (N)


          The scale factors for A.

SCOND

          SCOND is DOUBLE PRECISION


          Ratio of the smallest S(i) to the largest S(i).

AMAX

          AMAX is DOUBLE PRECISION


          Absolute value of largest matrix entry.

EQUED

          EQUED is CHARACTER*1


          Specifies whether or not equilibration was done.


          = 'N':  No equilibration.


          = 'Y':  Equilibration was done, i.e., A has been replaced by


                  diag(S) * A * diag(S).

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


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


  scaling is done.


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


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


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAQSP equilibrates a symmetric matrix A using the scaling factors


 in the vector S.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          symmetric matrix A is stored.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

N

          N is INTEGER


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

AP

          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)


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


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


          is stored in the array AP as follows:


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


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


          On exit, the equilibrated matrix:  diag(S) * A * diag(S), in


          the same storage format as A.

S

          S is DOUBLE PRECISION array, dimension (N)


          The scale factors for A.

SCOND

          SCOND is DOUBLE PRECISION


          Ratio of the smallest S(i) to the largest S(i).

AMAX

          AMAX is DOUBLE PRECISION


          Absolute value of largest matrix entry.

EQUED

          EQUED is CHARACTER*1


          Specifies whether or not equilibration was done.


          = 'N':  No equilibration.


          = 'Y':  Equilibration was done, i.e., A has been replaced by


                  diag(S) * A * diag(S).

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


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


  scaling is done.


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


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


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAQTR solves the real quasi-triangular system


              op(T)*p = scale*c,               if LREAL = .TRUE.


 or the complex quasi-triangular systems


            op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.


 in real arithmetic, where T is upper quasi-triangular.


 If LREAL = .FALSE., then the first diagonal block of T must be


 1 by 1, B is the specially structured matrix


                B = [ b(1) b(2) ... b(n) ]


                    [       w            ]


                    [           w        ]


                    [              .     ]


                    [                 w  ]


 op(A) = A or A**T, A**T denotes the transpose of


 matrix A.


 On input, X = [ c ].  On output, X = [ p ].


               [ d ]                  [ q ]


 This subroutine is designed for the condition number estimation


 in routine DTRSNA.

LTRAN

          LTRAN is LOGICAL


          On entry, LTRAN specifies the option of conjugate transpose:


             = .FALSE.,    op(T+i*B) = T+i*B,


             = .TRUE.,     op(T+i*B) = (T+i*B)**T.

LREAL

          LREAL is LOGICAL


          On entry, LREAL specifies the input matrix structure:


             = .FALSE.,    the input is complex


             = .TRUE.,     the input is real

N

          N is INTEGER


          On entry, N specifies the order of T+i*B. N >= 0.

T

          T is DOUBLE PRECISION array, dimension (LDT,N)


          On entry, T contains a matrix in Schur canonical form.


          If LREAL = .FALSE., then the first diagonal block of T mu


          be 1 by 1.

LDT

          LDT is INTEGER


          The leading dimension of the matrix T. LDT >= max(1,N).

B

          B is DOUBLE PRECISION array, dimension (N)


          On entry, B contains the elements to form the matrix


          B as described above.


          If LREAL = .TRUE., B is not referenced.

W

          W is DOUBLE PRECISION


          On entry, W is the diagonal element of the matrix B.


          If LREAL = .TRUE., W is not referenced.

SCALE

          SCALE is DOUBLE PRECISION


          On exit, SCALE is the scale factor.

X

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


          On entry, X contains the right hand side of the system.


          On exit, X is overwritten by the solution.

WORK

          WORK is DOUBLE PRECISION array, dimension (N)

INFO

          INFO is INTEGER


          On exit, INFO is set to


             0: successful exit.


               1: the some diagonal 1 by 1 block has been perturbed by


                  a small number SMIN to keep nonsingularity.


               2: the some diagonal 2 by 2 block has been perturbed by


                  a small number in DLALN2 to keep nonsingularity.


          NOTE: In the interests of speed, this routine does not


                check the inputs for errors.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAR1V computes the (scaled) r-th column of the inverse of


 the sumbmatrix in rows B1 through BN of the tridiagonal matrix


 L D L**T - sigma I. When sigma is close to an eigenvalue, the


 computed vector is an accurate eigenvector. Usually, r corresponds


 to the index where the eigenvector is largest in magnitude.


 The following steps accomplish this computation :


 (a) Stationary qd transform,  L D L**T - sigma I = L(+) D(+) L(+)**T,


 (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,


 (c) Computation of the diagonal elements of the inverse of


     L D L**T - sigma I by combining the above transforms, and choosing


     r as the index where the diagonal of the inverse is (one of the)


     largest in magnitude.


 (d) Computation of the (scaled) r-th column of the inverse using the


     twisted factorization obtained by combining the top part of the


     the stationary and the bottom part of the progressive transform.

N

          N is INTEGER


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

B1

          B1 is INTEGER


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

BN

          BN is INTEGER


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

LAMBDA

          LAMBDA is DOUBLE PRECISION


           The shift. In order to compute an accurate eigenvector,


           LAMBDA should be a good approximation to an eigenvalue


           of L D L**T.

L

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


           The (n-1) subdiagonal elements of the unit bidiagonal matrix


           L, in elements 1 to N-1.

D

          D is DOUBLE PRECISION array, dimension (N)


           The n diagonal elements of the diagonal matrix D.

LD

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


           The n-1 elements L(i)*D(i).

LLD

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


           The n-1 elements L(i)*L(i)*D(i).

PIVMIN

          PIVMIN is DOUBLE PRECISION


           The minimum pivot in the Sturm sequence.

GAPTOL

          GAPTOL is DOUBLE PRECISION


           Tolerance that indicates when eigenvector entries are negligible


           w.r.t. their contribution to the residual.

Z

          Z is DOUBLE PRECISION array, dimension (N)


           On input, all entries of Z must be set to 0.


           On output, Z contains the (scaled) r-th column of the


           inverse. The scaling is such that Z(R) equals 1.

WANTNC

          WANTNC is LOGICAL


           Specifies whether NEGCNT has to be computed.

NEGCNT

          NEGCNT is INTEGER


           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin


           in the  matrix factorization L D L**T, and NEGCNT = -1 otherwise.

ZTZ

          ZTZ is DOUBLE PRECISION


           The square of the 2-norm of Z.

MINGMA

          MINGMA is DOUBLE PRECISION


           The reciprocal of the largest (in magnitude) diagonal


           element of the inverse of L D L**T - sigma I.

R

          R is INTEGER


           The twist index for the twisted factorization used to


           compute Z.


           On input, 0 <= R <= N. If R is input as 0, R is set to


           the index where (L D L**T - sigma I)^{-1} is largest


           in magnitude. If 1 <= R <= N, R is unchanged.


           On output, R contains the twist index used to compute Z.


           Ideally, R designates the position of the maximum entry in the


           eigenvector.

ISUPPZ

          ISUPPZ is INTEGER array, dimension (2)


           The support of the vector in Z, i.e., the vector Z is


           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).

NRMINV

          NRMINV is DOUBLE PRECISION


           NRMINV = 1/SQRT( ZTZ )

RESID

          RESID is DOUBLE PRECISION


           The residual of the FP vector.


           RESID = ABS( MINGMA )/SQRT( ZTZ )

RQCORR

          RQCORR is DOUBLE PRECISION


           The Rayleigh Quotient correction to LAMBDA.


           RQCORR = MINGMA*TMP

WORK

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

 DLAR2V applies a vector of real plane rotations from both sides to


 a sequence of 2-by-2 real symmetric matrices, defined by the elements


 of the vectors x, y and z. For i = 1,2,...,n


    ( x(i)  z(i) ) := (  c(i)  s(i) ) ( x(i)  z(i) ) ( c(i) -s(i) )


    ( z(i)  y(i) )    ( -s(i)  c(i) ) ( z(i)  y(i) ) ( s(i)  c(i) )

N

          N is INTEGER


          The number of plane rotations to be applied.

X

          X is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCX)


          The vector x.

Y

          Y is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCX)


          The vector y.

Z

          Z is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCX)


          The vector z.

INCX

          INCX is INTEGER


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

C

          C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)


          The cosines of the plane rotations.

S

          S is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)


          The sines of the plane rotations.

INCC

          INCC is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARF applies a real elementary reflector H to a real m by n matrix


 C, from either the left or the right. H is represented in the form


       H = I - tau * v * v**T


 where tau is a real scalar and v is a real vector.


 If tau = 0, then H is taken to be the unit matrix.

SIDE

          SIDE is CHARACTER*1


          = 'L': form  H * C


          = 'R': form  C * H

M

          M is INTEGER


          The number of rows of the matrix C.

N

          N is INTEGER


          The number of columns of the matrix C.

V

          V is DOUBLE PRECISION array, dimension


                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'


                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'


          The vector v in the representation of H. V is not used if


          TAU = 0.

INCV

          INCV is INTEGER


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

TAU

          TAU is DOUBLE PRECISION


          The value tau in the representation of H.

C

          C is DOUBLE PRECISION array, dimension (LDC,N)


          On entry, the m by n matrix C.


          On exit, C is overwritten by the matrix H * C if SIDE = 'L',


          or C * H if SIDE = 'R'.

LDC

          LDC is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension


                         (N) if SIDE = 'L'


                      or (M) if SIDE = 'R'

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARFB applies a real block reflector H or its transpose H**T to a


 real m by n matrix C, from either the left or the right.

SIDE

          SIDE is CHARACTER*1


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


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

TRANS

          TRANS is CHARACTER*1


          = 'N': apply H (No transpose)


          = 'T': apply H**T (Transpose)

DIRECT

          DIRECT is CHARACTER*1


          Indicates how H is formed from a product of elementary


          reflectors


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


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

STOREV

          STOREV is CHARACTER*1


          Indicates how the vectors which define the elementary


          reflectors are stored:


          = 'C': Columnwise


          = 'R': Rowwise

M

          M is INTEGER


          The number of rows of the matrix C.

N

          N is INTEGER


          The number of columns of the matrix C.

K

          K is INTEGER


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


          reflectors whose product defines the block reflector).


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


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

V

          V is DOUBLE PRECISION array, dimension


                                (LDV,K) if STOREV = 'C'


                                (LDV,M) if STOREV = 'R' and SIDE = 'L'


                                (LDV,N) if STOREV = 'R' and SIDE = 'R'


          The matrix V. See Further Details.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);


          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);


          if STOREV = 'R', LDV >= K.

T

          T is DOUBLE PRECISION array, dimension (LDT,K)


          The triangular k by k matrix T in the representation of the


          block reflector.

LDT

          LDT is INTEGER


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

C

          C is DOUBLE PRECISION array, dimension (LDC,N)


          On entry, the m by n matrix C.


          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.

LDC

          LDC is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)

LDWORK

          LDWORK is INTEGER


          The leading dimension of the array WORK.


          If SIDE = 'L', LDWORK >= max(1,N);


          if SIDE = 'R', LDWORK >= max(1,M).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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


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


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


  array is not used.


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


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


                   ( v1  1    )                     (     1 v2 v2 v2 )


                   ( v1 v2  1 )                     (        1 v3 v3 )


                   ( v1 v2 v3 )


                   ( v1 v2 v3 )


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


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


                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )


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


                   (     1 v3 )


                   (        1 )

 DLARFB_GETT applies a real Householder block reflector H from the


 left to a real (K+M)-by-N  'triangular-pentagonal' matrix


 composed of two block matrices: an upper trapezoidal K-by-N matrix A


 stored in the array A, and a rectangular M-by-(N-K) matrix B, stored


 in the array B. The block reflector H is stored in a compact


 WY-representation, where the elementary reflectors are in the


 arrays A, B and T. See Further Details section.

IDENT

          IDENT is CHARACTER*1


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


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


             the input matrix A,


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


             not stored.


          See Further Details section.

M

          M is INTEGER


          The number of rows of the matrix B.


          M >= 0.

N

          N is INTEGER


          The number of columns of the matrices A and B.


          N >= 0.

K

          K is INTEGER


          The number or rows of the matrix A.


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


          elementary reflectors whose product defines the block


          reflector. 0 <= K <= N.

T

          T is DOUBLE PRECISION array, dimension (LDT,K)


          The upper-triangular K-by-K matrix T in the representation


          of the block reflector.

LDT

          LDT is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry:


           a) In the K-by-N upper-trapezoidal part A: input matrix A.


           b) In the columns below the diagonal: columns of V1


              (ones are not stored on the diagonal).


          On exit:


            A is overwritten by rectangular K-by-N product H*A.


          See Further Details section.

LDA

          LDB is INTEGER


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

B

          B is DOUBLE PRECISION array, dimension (LDB,N)


          On entry:


            a) In the M-by-(N-K) right block: input matrix B.


            b) In the M-by-N left block: columns of V2.


          On exit:


            B is overwritten by rectangular M-by-N product H*B.


          See Further Details section.

LDB

          LDB is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array,


          dimension (LDWORK,max(K,N-K))

LDWORK

          LDWORK is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 November 2020, Igor Kozachenko,


                Computer Science Division,


                University of California, Berkeley

    (1) Description of the Algebraic Operation.


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


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


    A = ( A1, A2 ).


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


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


    B = ( B1, B2 ).


    Perform the operation:


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


       ( B_out )        ( B_in )                          ( B_in )


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


                          ( V2 )                            ( B_in )


     On input:


    a) ( A_in )  consists of two block columns:


       ( B_in )


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


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


       where the column blocks are:


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


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


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


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


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


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


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


    b) V = ( V1 )


           ( V2 )


       where:


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


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


          stored in the lower-triangular part of the array


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


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


                 (because on input B1_in is a rectangular zero


                  matrix that is not stored and the space is


                  used to store V2).


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


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


    On output:


    a) ( A_out ) consists of two  block columns:


       ( B_out )


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


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


       where the column blocks are:


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


                   upper-triangular matrix, if V1 is an


                   identity matrix. AiOut is stored in


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


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


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


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


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


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


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


    The operation above can be represented as the same operation


    on each block column:


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


       ( B1_out )        (     0 )                          (     0 )


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


       ( B2_out )        ( B2_in )                          ( B2_in )


    If IDENT != 'I':


       The computation for column block 1:


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


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


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


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


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


    If IDENT == 'I':


       The operation for column block 1:


       A1_out: = A1_in - V1*T**A1_in


       B1_out: = - V2*T**A1_in


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


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


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


    (2) Description of the Algorithmic Computation.


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


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


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


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


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


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


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


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


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


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


    the first step and W1 from the second step.


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


    more computations than in the Case (B).


    if( IDENT != 'I' ) then


     if ( N > K ) then


       (First Step - column block 2)


       col2_(1) W2: = A2


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


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


       col2_(4) W2: = T * W2


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


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


       col2_(7) A2: = A2 - W2


     else


       (Second Step - column block 1)


       col1_(1) W1: = A1


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


       col1_(3) W1: = T * W1


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


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


       col1_(6) square A1: = A1 - W1


     end if


    end if


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


    less computations than in the Case (A)


    if( IDENT == 'I' ) then


     if ( N > K ) then


       (First Step - column block 2)


       col2_(1) W2: = A2


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


       col2_(4) W2: = T * W2


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


       col2_(7) A2: = A2 - W2


     else


       (Second Step - column block 1)


       col1_(1) W1: = A1


       col1_(3) W1: = T * W1


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


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


     end if


    end if


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


    algorithm:


    if ( N > K ) then


      (First Step - column block 2)


      col2_(1)  W2: = A2


      if( IDENT != 'I' ) then


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


                      = (unit_lower_tr_of_(A1)**T) * W2


      end if


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


      col2_(4)  W2: = T * W2


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


      if( IDENT != 'I' ) then


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


      end if


      col2_(7) A2: = A2 - W2


    else


    (Second Step - column block 1)


      col1_(1) W1: = A1


      if( IDENT != 'I' ) then


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


                    = (unit_lower_tr_of_(A1)**T) * W1


      end if


      col1_(3) W1: = T * W1


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


      if( IDENT != 'I' ) then


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


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


      end if


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


    end if

 DLARFG generates a real elementary reflector H of order n, such


 that


       H * ( alpha ) = ( beta ),   H**T * H = I.


           (   x   )   (   0  )


 where alpha and beta are scalars, and x is an (n-1)-element real


 vector. H is represented in the form


       H = I - tau * ( 1 ) * ( 1 v**T ) ,


                     ( v )


 where tau is a real scalar and v is a real (n-1)-element


 vector.


 If the elements of x are all zero, then tau = 0 and H is taken to be


 the unit matrix.


 Otherwise  1 <= tau <= 2.

N

          N is INTEGER


          The order of the elementary reflector.

ALPHA

          ALPHA is DOUBLE PRECISION


          On entry, the value alpha.


          On exit, it is overwritten with the value beta.

X

          X is DOUBLE PRECISION array, dimension


                         (1+(N-2)*abs(INCX))


          On entry, the vector x.


          On exit, it is overwritten with the vector v.

INCX

          INCX is INTEGER


          The increment between elements of X. INCX > 0.

TAU

          TAU is DOUBLE PRECISION


          The value tau.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARFGP generates a real elementary reflector H of order n, such


 that


       H * ( alpha ) = ( beta ),   H**T * H = I.


           (   x   )   (   0  )


 where alpha and beta are scalars, beta is non-negative, and x is


 an (n-1)-element real vector.  H is represented in the form


       H = I - tau * ( 1 ) * ( 1 v**T ) ,


                     ( v )


 where tau is a real scalar and v is a real (n-1)-element


 vector.


 If the elements of x are all zero, then tau = 0 and H is taken to be


 the unit matrix.

N

          N is INTEGER


          The order of the elementary reflector.

ALPHA

          ALPHA is DOUBLE PRECISION


          On entry, the value alpha.


          On exit, it is overwritten with the value beta.

X

          X is DOUBLE PRECISION array, dimension


                         (1+(N-2)*abs(INCX))


          On entry, the vector x.


          On exit, it is overwritten with the vector v.

INCX

          INCX is INTEGER


          The increment between elements of X. INCX > 0.

TAU

          TAU is DOUBLE PRECISION


          The value tau.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARFT forms the triangular factor T of a real block reflector H


 of order n, which is defined as a product of k elementary reflectors.


 If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;


 If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.


 If STOREV = 'C', the vector which defines the elementary reflector


 H(i) is stored in the i-th column of the array V, and


    H  =  I - V * T * V**T


 If STOREV = 'R', the vector which defines the elementary reflector


 H(i) is stored in the i-th row of the array V, and


    H  =  I - V**T * T * V

DIRECT

          DIRECT is CHARACTER*1


          Specifies the order in which the elementary reflectors are


          multiplied to form the block reflector:


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


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

STOREV

          STOREV is CHARACTER*1


          Specifies how the vectors which define the elementary


          reflectors are stored (see also Further Details):


          = 'C': columnwise


          = 'R': rowwise

N

          N is INTEGER


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

K

          K is INTEGER


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


          elementary reflectors). K >= 1.

V

          V is DOUBLE PRECISION array, dimension


                               (LDV,K) if STOREV = 'C'


                               (LDV,N) if STOREV = 'R'


          The matrix V. See further details.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.

TAU

          TAU is DOUBLE PRECISION array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i).

T

          T is DOUBLE PRECISION array, dimension (LDT,K)


          The k by k triangular factor T of the block reflector.


          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is


          lower triangular. The rest of the array is not used.

LDT

          LDT is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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


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


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


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


                   ( v1  1    )                     (     1 v2 v2 v2 )


                   ( v1 v2  1 )                     (        1 v3 v3 )


                   ( v1 v2 v3 )


                   ( v1 v2 v3 )


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


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


                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )


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


                   (     1 v3 )


                   (        1 )

 DLARFX applies a real elementary reflector H to a real m by n


 matrix C, from either the left or the right. H is represented in the


 form


       H = I - tau * v * v**T


 where tau is a real scalar and v is a real vector.


 If tau = 0, then H is taken to be the unit matrix


 This version uses inline code if H has order < 11.

SIDE

          SIDE is CHARACTER*1


          = 'L': form  H * C


          = 'R': form  C * H

M

          M is INTEGER


          The number of rows of the matrix C.

N

          N is INTEGER


          The number of columns of the matrix C.

V

          V is DOUBLE PRECISION array, dimension (M) if SIDE = 'L'


                                     or (N) if SIDE = 'R'


          The vector v in the representation of H.

TAU

          TAU is DOUBLE PRECISION


          The value tau in the representation of H.

C

          C is DOUBLE PRECISION array, dimension (LDC,N)


          On entry, the m by n matrix C.


          On exit, C is overwritten by the matrix H * C if SIDE = 'L',


          or C * H if SIDE = 'R'.

LDC

          LDC is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension


                      (N) if SIDE = 'L'


                      or (M) if SIDE = 'R'


          WORK is not referenced if H has order < 11.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARFY applies an elementary reflector, or Householder matrix, H,


 to an n x n symmetric matrix C, from both the left and the right.


 H is represented in the form


    H = I - tau * v * v'


 where  tau  is a scalar and  v  is a vector.


 If  tau  is  zero, then  H  is taken to be the unit matrix.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          symmetric matrix C is stored.


          = 'U':  Upper triangle


          = 'L':  Lower triangle

N

          N is INTEGER


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

V

          V is DOUBLE PRECISION array, dimension


                  (1 + (N-1)*abs(INCV))


          The vector v as described above.

INCV

          INCV is INTEGER


          The increment between successive elements of v.  INCV must


          not be zero.

TAU

          TAU is DOUBLE PRECISION


          The value tau as described above.

C

          C is DOUBLE PRECISION array, dimension (LDC, N)


          On entry, the matrix C.


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

LDC

          LDC is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (N)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARGV generates a vector of real plane rotations, determined by


 elements of the real vectors x and y. For i = 1,2,...,n


    (  c(i)  s(i) ) ( x(i) ) = ( a(i) )


    ( -s(i)  c(i) ) ( y(i) ) = (   0  )

N

          N is INTEGER


          The number of plane rotations to be generated.

X

          X is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCX)


          On entry, the vector x.


          On exit, x(i) is overwritten by a(i), for i = 1,...,n.

INCX

          INCX is INTEGER


          The increment between elements of X. INCX > 0.

Y

          Y is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCY)


          On entry, the vector y.


          On exit, the sines of the plane rotations.

INCY

          INCY is INTEGER


          The increment between elements of Y. INCY > 0.

C

          C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)


          The cosines of the plane rotations.

INCC

          INCC is INTEGER


          The increment between elements of C. INCC > 0.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLARRV computes the eigenvectors of the tridiagonal matrix


 T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T.


 The input eigenvalues should have been computed by DLARRE.

N

          N is INTEGER


          The order of the matrix.  N >= 0.

VL

          VL is DOUBLE PRECISION


          Lower bound of the interval that contains the desired


          eigenvalues. VL < VU. Needed to compute gaps on the left or right


          end of the extremal eigenvalues in the desired RANGE.

VU

          VU is DOUBLE PRECISION


          Upper bound of the interval that contains the desired


          eigenvalues. VL < VU. 


          Note: VU is currently not used by this implementation of DLARRV, VU is


          passed to DLARRV because it could be used compute gaps on the right end


          of the extremal eigenvalues. However, with not much initial accuracy in


          LAMBDA and VU, the formula can lead to an overestimation of the right gap


          and thus to inadequately early RQI 'convergence'. This is currently


          prevented this by forcing a small right gap. And so it turns out that VU


          is currently not used by this implementation of DLARRV.

D

          D is DOUBLE PRECISION array, dimension (N)


          On entry, the N diagonal elements of the diagonal matrix D.


          On exit, D may be overwritten.

L

          L is DOUBLE PRECISION array, dimension (N)


          On entry, the (N-1) subdiagonal elements of the unit


          bidiagonal matrix L are in elements 1 to N-1 of L


          (if the matrix is not split.) At the end of each block


          is stored the corresponding shift as given by DLARRE.


          On exit, L is overwritten.

PIVMIN

          PIVMIN is DOUBLE PRECISION


          The minimum pivot allowed in the Sturm sequence.

ISPLIT

          ISPLIT is INTEGER array, dimension (N)


          The splitting points, at which T breaks up into blocks.


          The first block consists of rows/columns 1 to


          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1


          through ISPLIT( 2 ), etc.

M

          M is INTEGER


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

DOL

          DOL is INTEGER

DOU

          DOU is INTEGER


          If the user wants to compute only selected eigenvectors from all


          the eigenvalues supplied, he can specify an index range DOL:DOU.


          Or else the setting DOL=1, DOU=M should be applied.


          Note that DOL and DOU refer to the order in which the eigenvalues


          are stored in W.


          If the user wants to compute only selected eigenpairs, then


          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the


          computed eigenvectors. All other columns of Z are set to zero.

MINRGP

          MINRGP is DOUBLE PRECISION

RTOL1

          RTOL1 is DOUBLE PRECISION

RTOL2

          RTOL2 is DOUBLE PRECISION


           Parameters for bisection.


           An interval [LEFT,RIGHT] has converged if


           RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )

W

          W is DOUBLE PRECISION array, dimension (N)


          The first M elements of W contain the APPROXIMATE eigenvalues for


          which eigenvectors are to be computed.  The eigenvalues


          should be grouped by split-off block and ordered from


          smallest to largest within the block ( The output array


          W from DLARRE is expected here ). Furthermore, they are with


          respect to the shift of the corresponding root representation


          for their block. On exit, W holds the eigenvalues of the


          UNshifted matrix.

WERR

          WERR is DOUBLE PRECISION array, dimension (N)


          The first M elements contain the semiwidth of the uncertainty


          interval of the corresponding eigenvalue in W

WGAP

          WGAP is DOUBLE PRECISION array, dimension (N)


          The separation from the right neighbor eigenvalue in W.

IBLOCK

          IBLOCK is INTEGER array, dimension (N)


          The indices of the blocks (submatrices) associated with the


          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue


          W(i) belongs to the first block from the top, =2 if W(i)


          belongs to the second block, etc.

INDEXW

          INDEXW is INTEGER array, dimension (N)


          The indices of the eigenvalues within each block (submatrix);


          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the


          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block.

GERS

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


          The N Gerschgorin intervals (the i-th Gerschgorin interval


          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should


          be computed from the original UNshifted matrix.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )


          If INFO = 0, the first M columns of Z contain the


          orthonormal eigenvectors of the matrix T


          corresponding to the input eigenvalues, with the i-th


          column of Z holding the eigenvector associated with W(i).


          Note: the user must ensure that at least max(1,M) columns are


          supplied in the array Z.

LDZ

          LDZ is INTEGER


          The leading dimension of the array Z.  LDZ >= 1, and if


          JOBZ = 'V', LDZ >= max(1,N).

ISUPPZ

          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )


          The support of the eigenvectors in Z, i.e., the indices


          indicating the nonzero elements in Z. The I-th eigenvector


          is nonzero only in elements ISUPPZ( 2*I-1 ) through


          ISUPPZ( 2*I ).

WORK

          WORK is DOUBLE PRECISION array, dimension (12*N)

IWORK

          IWORK is INTEGER array, dimension (7*N)

INFO

          INFO is INTEGER


          = 0:  successful exit


          > 0:  A problem occurred in DLARRV.


          < 0:  One of the called subroutines signaled an internal problem.


                Needs inspection of the corresponding parameter IINFO


                for further information.


          =-1:  Problem in DLARRB when refining a child's eigenvalues.


          =-2:  Problem in DLARRF when computing the RRR of a child.


                When a child is inside a tight cluster, it can be difficult


                to find an RRR. A partial remedy from the user's point of


                view is to make the parameter MINRGP smaller and recompile.


                However, as the orthogonality of the computed vectors is


                proportional to 1/MINRGP, the user should be aware that


                he might be trading in precision when he decreases MINRGP.


          =-3:  Problem in DLARRB when refining a single eigenvalue


                after the Rayleigh correction was rejected.


          = 5:  The Rayleigh Quotient Iteration failed to converge to


                full accuracy in MAXITR steps.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

 DLARTV applies a vector of real plane rotations to elements of the


 real vectors x and y. For i = 1,2,...,n


    ( x(i) ) := (  c(i)  s(i) ) ( x(i) )


    ( y(i) )    ( -s(i)  c(i) ) ( y(i) )

N

          N is INTEGER


          The number of plane rotations to be applied.

X

          X is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCX)


          The vector x.

INCX

          INCX is INTEGER


          The increment between elements of X. INCX > 0.

Y

          Y is DOUBLE PRECISION array,


                         dimension (1+(N-1)*INCY)


          The vector y.

INCY

          INCY is INTEGER


          The increment between elements of Y. INCY > 0.

C

          C is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)


          The cosines of the plane rotations.

S

          S is DOUBLE PRECISION array, dimension (1+(N-1)*INCC)


          The sines of the plane rotations.

INCC

          INCC is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLASWP performs a series of row interchanges on the matrix A.


 One row interchange is initiated for each of rows K1 through K2 of A.

N

          N is INTEGER


          The number of columns of the matrix A.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the matrix of column dimension N to which the row


          interchanges will be applied.


          On exit, the permuted matrix.

LDA

          LDA is INTEGER


          The leading dimension of the array A.

K1

          K1 is INTEGER


          The first element of IPIV for which a row interchange will


          be done.

K2

          K2 is INTEGER


          (K2-K1+1) is the number of elements of IPIV for which a row


          interchange will be done.

IPIV

          IPIV is INTEGER array, dimension (K1+(K2-K1)*abs(INCX))


          The vector of pivot indices. Only the elements in positions


          K1 through K1+(K2-K1)*abs(INCX) of IPIV are accessed.


          IPIV(K1+(K-K1)*abs(INCX)) = L implies rows K and L are to be


          interchanged.

INCX

          INCX is INTEGER


          The increment between successive values of IPIV. If INCX


          is negative, the pivots are applied in reverse order.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Modified by


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

 DLAT2S converts a DOUBLE PRECISION triangular matrix, SA, to a SINGLE


 PRECISION triangular matrix, A.


 RMAX is the overflow for the SINGLE PRECISION arithmetic


 DLAS2S checks that all the entries of A are between -RMAX and


 RMAX. If not the conversion is aborted and a flag is raised.


 This is an auxiliary routine so there is no argument checking.

UPLO

          UPLO is CHARACTER*1


          = 'U':  A is upper triangular;


          = 'L':  A is lower triangular.

N

          N is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the N-by-N triangular coefficient matrix A.

LDA

          LDA is INTEGER


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

SA

          SA is REAL array, dimension (LDSA,N)


          Only the UPLO part of SA is referenced.  On exit, if INFO=0,


          the N-by-N coefficient matrix SA; if INFO>0, the content of


          the UPLO part of SA is unspecified.

LDSA

          LDSA is INTEGER


          The leading dimension of the array SA.  LDSA >= max(1,M).

INFO

          INFO is INTEGER


          = 0:  successful exit.


          = 1:  an entry of the matrix A is greater than the SINGLE


                PRECISION overflow threshold, in this case, the content


                of the UPLO part of SA in exit is unspecified.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLATBS solves one of the triangular systems


    A *x = s*b  or  A**T*x = s*b


 with scaling to prevent overflow, where A is an upper or lower


 triangular band matrix.  Here A**T denotes the transpose of A, x and b


 are n-element vectors, and s is a scaling factor, usually less than


 or equal to 1, chosen so that the components of x will be less than


 the overflow threshold.  If the unscaled problem will not cause


 overflow, the Level 2 BLAS routine DTBSV is called.  If the matrix A


 is singular (A(j,j) = 0 for some j), then s is set to 0 and a


 non-trivial solution to A*x = 0 is returned.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


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

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

NORMIN

          NORMIN is CHARACTER*1


          Specifies whether CNORM has been set or not.


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


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


                  be computed and stored in CNORM.

N

          N is INTEGER


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

KD

          KD is INTEGER


          The number of subdiagonals or superdiagonals in the


          triangular matrix A.  KD >= 0.

AB

          AB is DOUBLE PRECISION array, dimension (LDAB,N)


          The upper or lower triangular band matrix A, stored in the


          first KD+1 rows of the array. The j-th column of A is stored


          in the j-th column of the array AB as follows:


          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;


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

LDAB

          LDAB is INTEGER


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

X

          X is DOUBLE PRECISION array, dimension (N)


          On entry, the right hand side b of the triangular system.


          On exit, X is overwritten by the solution vector x.

SCALE

          SCALE is DOUBLE PRECISION


          The scaling factor s for the triangular system


             A * x = s*b  or  A**T* x = s*b.


          If SCALE = 0, the matrix A is singular or badly scaled, and


          the vector x is an exact or approximate solution to A*x = 0.

CNORM

          CNORM is DOUBLE PRECISION array, dimension (N)


          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)


          contains the norm of the off-diagonal part of the j-th column


          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal


          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)


          must be greater than or equal to the 1-norm.


          If NORMIN = 'N', CNORM is an output argument and CNORM(j)


          returns the 1-norm of the offdiagonal part of the j-th column


          of A.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  A rough bound on x is computed; if that is less than overflow, DTBSV


  is called, otherwise, specific code is used which checks for possible


  overflow or divide-by-zero at every operation.


  A columnwise scheme is used for solving A*x = b.  The basic algorithm


  if A is lower triangular is


       x[1:n] := b[1:n]


       for j = 1, ..., n


            x(j) := x(j) / A(j,j)


            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]


       end


  Define bounds on the components of x after j iterations of the loop:


     M(j) = bound on x[1:j]


     G(j) = bound on x[j+1:n]


  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.


  Then for iteration j+1 we have


     M(j+1) <= G(j) / | A(j+1,j+1) |


     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |


            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )


  where CNORM(j+1) is greater than or equal to the infinity-norm of


  column j+1 of A, not counting the diagonal.  Hence


     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )


                  1<=i<=j


  and


     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )


                                   1<=i< j


  Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the


  reciprocal of the largest M(j), j=1,..,n, is larger than


  max(underflow, 1/overflow).


  The bound on x(j) is also used to determine when a step in the


  columnwise method can be performed without fear of overflow.  If


  the computed bound is greater than a large constant, x is scaled to


  prevent overflow, but if the bound overflows, x is set to 0, x(j) to


  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.


  Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic


  algorithm for A upper triangular is


       for j = 1, ..., n


            x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)


       end


  We simultaneously compute two bounds


       G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j


       M(j) = bound on x(i), 1<=i<=j


  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we


  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.


  Then the bound on x(j) is


       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |


            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )


                      1<=i<=j


  and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater


  than max(underflow, 1/overflow).

 DLATDF uses the LU factorization of the n-by-n matrix Z computed by


 DGETC2 and computes a contribution to the reciprocal Dif-estimate


 by solving Z * x = b for x, and choosing the r.h.s. b such that


 the norm of x is as large as possible. On entry RHS = b holds the


 contribution from earlier solved sub-systems, and on return RHS = x.


 The factorization of Z returned by DGETC2 has the form Z = P*L*U*Q,


 where P and Q are permutation matrices. L is lower triangular with


 unit diagonal elements and U is upper triangular.

IJOB

          IJOB is INTEGER


          IJOB = 2: First compute an approximative null-vector e


              of Z using DGECON, e is normalized and solve for


              Zx = +-e - f with the sign giving the greater value


              of 2-norm(x). About 5 times as expensive as Default.


          IJOB .ne. 2: Local look ahead strategy where all entries of


              the r.h.s. b is chosen as either +1 or -1 (Default).

N

          N is INTEGER


          The number of columns of the matrix Z.

Z

          Z is DOUBLE PRECISION array, dimension (LDZ, N)


          On entry, the LU part of the factorization of the n-by-n


          matrix Z computed by DGETC2:  Z = P * L * U * Q

LDZ

          LDZ is INTEGER


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

RHS

          RHS is DOUBLE PRECISION array, dimension (N)


          On entry, RHS contains contributions from other subsystems.


          On exit, RHS contains the solution of the subsystem with


          entries according to the value of IJOB (see above).

RDSUM

          RDSUM is DOUBLE PRECISION


          On entry, the sum of squares of computed contributions to


          the Dif-estimate under computation by DTGSYL, where the


          scaling factor RDSCAL (see below) has been factored out.


          On exit, the corresponding sum of squares updated with the


          contributions from the current sub-system.


          If TRANS = 'T' RDSUM is not touched.


          NOTE: RDSUM only makes sense when DTGSY2 is called by STGSYL.

RDSCAL

          RDSCAL is DOUBLE PRECISION


          On entry, scaling factor used to prevent overflow in RDSUM.


          On exit, RDSCAL is updated w.r.t. the current contributions


          in RDSUM.


          If TRANS = 'T', RDSCAL is not touched.


          NOTE: RDSCAL only makes sense when DTGSY2 is called by


                DTGSYL.

IPIV

          IPIV is INTEGER array, dimension (N).


          The pivot indices; for 1 <= i <= N, row i of the


          matrix has been interchanged with row IPIV(i).

JPIV

          JPIV is INTEGER array, dimension (N).


          The pivot indices; for 1 <= j <= N, column j of the


          matrix has been interchanged with column JPIV(j).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

  [1] Bo Kagstrom and Lars Westin,


      Generalized Schur Methods with Condition Estimators for


      Solving the Generalized Sylvester Equation, IEEE Transactions


      on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.


  [2] Peter Poromaa,


      On Efficient and Robust Estimators for the Separation


      between two Regular Matrix Pairs with Applications in


      Condition Estimation. Report IMINF-95.05, Departement of


      Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.

 DLATPS solves one of the triangular systems


    A *x = s*b  or  A**T*x = s*b


 with scaling to prevent overflow, where A is an upper or lower


 triangular matrix stored in packed form.  Here A**T denotes the


 transpose of A, x and b are n-element vectors, and s is a scaling


 factor, usually less than or equal to 1, chosen so that the


 components of x will be less than the overflow threshold.  If the


 unscaled problem will not cause overflow, the Level 2 BLAS routine


 DTPSV is called. If the matrix A is singular (A(j,j) = 0 for some j),


 then s is set to 0 and a non-trivial solution to A*x = 0 is returned.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


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

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

NORMIN

          NORMIN is CHARACTER*1


          Specifies whether CNORM has been set or not.


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


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


                  be computed and stored in CNORM.

N

          N is INTEGER


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

AP

          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)


          The upper or lower triangular matrix A, packed columnwise in


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


          AP as follows:


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


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

X

          X is DOUBLE PRECISION array, dimension (N)


          On entry, the right hand side b of the triangular system.


          On exit, X is overwritten by the solution vector x.

SCALE

          SCALE is DOUBLE PRECISION


          The scaling factor s for the triangular system


             A * x = s*b  or  A**T* x = s*b.


          If SCALE = 0, the matrix A is singular or badly scaled, and


          the vector x is an exact or approximate solution to A*x = 0.

CNORM

          CNORM is DOUBLE PRECISION array, dimension (N)


          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)


          contains the norm of the off-diagonal part of the j-th column


          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal


          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)


          must be greater than or equal to the 1-norm.


          If NORMIN = 'N', CNORM is an output argument and CNORM(j)


          returns the 1-norm of the offdiagonal part of the j-th column


          of A.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  A rough bound on x is computed; if that is less than overflow, DTPSV


  is called, otherwise, specific code is used which checks for possible


  overflow or divide-by-zero at every operation.


  A columnwise scheme is used for solving A*x = b.  The basic algorithm


  if A is lower triangular is


       x[1:n] := b[1:n]


       for j = 1, ..., n


            x(j) := x(j) / A(j,j)


            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]


       end


  Define bounds on the components of x after j iterations of the loop:


     M(j) = bound on x[1:j]


     G(j) = bound on x[j+1:n]


  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.


  Then for iteration j+1 we have


     M(j+1) <= G(j) / | A(j+1,j+1) |


     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |


            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )


  where CNORM(j+1) is greater than or equal to the infinity-norm of


  column j+1 of A, not counting the diagonal.  Hence


     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )


                  1<=i<=j


  and


     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )


                                   1<=i< j


  Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTPSV if the


  reciprocal of the largest M(j), j=1,..,n, is larger than


  max(underflow, 1/overflow).


  The bound on x(j) is also used to determine when a step in the


  columnwise method can be performed without fear of overflow.  If


  the computed bound is greater than a large constant, x is scaled to


  prevent overflow, but if the bound overflows, x is set to 0, x(j) to


  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.


  Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic


  algorithm for A upper triangular is


       for j = 1, ..., n


            x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)


       end


  We simultaneously compute two bounds


       G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j


       M(j) = bound on x(i), 1<=i<=j


  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we


  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.


  Then the bound on x(j) is


       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |


            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )


                      1<=i<=j


  and we can safely call DTPSV if 1/M(n) and 1/G(n) are both greater


  than max(underflow, 1/overflow).

 DLATRD reduces NB rows and columns of a real symmetric matrix A to


 symmetric tridiagonal form by an orthogonal similarity


 transformation Q**T * A * Q, and returns the matrices V and W which are


 needed to apply the transformation to the unreduced part of A.


 If UPLO = 'U', DLATRD reduces the last NB rows and columns of a


 matrix, of which the upper triangle is supplied;


 if UPLO = 'L', DLATRD reduces the first NB rows and columns of a


 matrix, of which the lower triangle is supplied.


 This is an auxiliary routine called by DSYTRD.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          symmetric matrix A is stored:


          = 'U': Upper triangular


          = 'L': Lower triangular

N

          N is INTEGER


          The order of the matrix A.

NB

          NB is INTEGER


          The number of rows and columns to be reduced.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the symmetric matrix A.  If UPLO = 'U', the leading


          n-by-n upper triangular part of A contains the upper


          triangular part of the matrix A, and the strictly lower


          triangular part of A is not referenced.  If UPLO = 'L', the


          leading n-by-n lower triangular part of A contains the lower


          triangular part of the matrix A, and the strictly upper


          triangular part of A is not referenced.


          On exit:


          if UPLO = 'U', the last NB columns have been reduced to


            tridiagonal form, with the diagonal elements overwriting


            the diagonal elements of A; the elements above the diagonal


            with the array TAU, represent the orthogonal matrix Q as a


            product of elementary reflectors;


          if UPLO = 'L', the first NB columns have been reduced to


            tridiagonal form, with the diagonal elements overwriting


            the diagonal elements of A; the elements below the diagonal


            with the array TAU, represent the  orthogonal matrix Q as a


            product of elementary reflectors.


          See Further Details.

LDA

          LDA is INTEGER


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

E

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


          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal


          elements of the last NB columns of the reduced matrix;


          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of


          the first NB columns of the reduced matrix.

TAU

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


          The scalar factors of the elementary reflectors, stored in


          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.


          See Further Details.

W

          W is DOUBLE PRECISION array, dimension (LDW,NB)


          The n-by-nb matrix W required to update the unreduced part


          of A.

LDW

          LDW is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


  reflectors


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


  Each H(i) has the form


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


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


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


  and tau in TAU(i-1).


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


  reflectors


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


  Each H(i) has the form


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


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


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


  and tau in TAU(i).


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


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


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


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


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


  with n = 5 and nb = 2:


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


    (  a   a   a   v4  v5 )              (  d                  )


    (      a   a   v4  v5 )              (  1   d              )


    (          a   1   v5 )              (  v1  1   a          )


    (              d   1  )              (  v1  v2  a   a      )


    (                  d  )              (  v1  v2  a   a   a  )


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


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


  an element of the vector defining H(i).

 DLATRS solves one of the triangular systems


    A *x = s*b  or  A**T *x = s*b


 with scaling to prevent overflow.  Here A is an upper or lower


 triangular matrix, A**T denotes the transpose of A, x and b are


 n-element vectors, and s is a scaling factor, usually less than


 or equal to 1, chosen so that the components of x will be less than


 the overflow threshold.  If the unscaled problem will not cause


 overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A


 is singular (A(j,j) = 0 for some j), then s is set to 0 and a


 non-trivial solution to A*x = 0 is returned.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


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

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

NORMIN

          NORMIN is CHARACTER*1


          Specifies whether CNORM has been set or not.


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


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


                  be computed and stored in CNORM.

N

          N is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          The triangular matrix A.  If UPLO = 'U', the leading n by n


          upper triangular part of the array A contains the upper


          triangular matrix, and the strictly lower triangular part of


          A is not referenced.  If UPLO = 'L', the leading n by n lower


          triangular part of the array A contains the lower triangular


          matrix, and the strictly upper triangular part of A is not


          referenced.  If DIAG = 'U', the diagonal elements of A are


          also not referenced and are assumed to be 1.

LDA

          LDA is INTEGER


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

X

          X is DOUBLE PRECISION array, dimension (N)


          On entry, the right hand side b of the triangular system.


          On exit, X is overwritten by the solution vector x.

SCALE

          SCALE is DOUBLE PRECISION


          The scaling factor s for the triangular system


             A * x = s*b  or  A**T* x = s*b.


          If SCALE = 0, the matrix A is singular or badly scaled, and


          the vector x is an exact or approximate solution to A*x = 0.

CNORM

          CNORM is DOUBLE PRECISION array, dimension (N)


          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)


          contains the norm of the off-diagonal part of the j-th column


          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal


          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)


          must be greater than or equal to the 1-norm.


          If NORMIN = 'N', CNORM is an output argument and CNORM(j)


          returns the 1-norm of the offdiagonal part of the j-th column


          of A.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  A rough bound on x is computed; if that is less than overflow, DTRSV


  is called, otherwise, specific code is used which checks for possible


  overflow or divide-by-zero at every operation.


  A columnwise scheme is used for solving A*x = b.  The basic algorithm


  if A is lower triangular is


       x[1:n] := b[1:n]


       for j = 1, ..., n


            x(j) := x(j) / A(j,j)


            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]


       end


  Define bounds on the components of x after j iterations of the loop:


     M(j) = bound on x[1:j]


     G(j) = bound on x[j+1:n]


  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.


  Then for iteration j+1 we have


     M(j+1) <= G(j) / | A(j+1,j+1) |


     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |


            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )


  where CNORM(j+1) is greater than or equal to the infinity-norm of


  column j+1 of A, not counting the diagonal.  Hence


     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )


                  1<=i<=j


  and


     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )


                                   1<=i< j


  Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the


  reciprocal of the largest M(j), j=1,..,n, is larger than


  max(underflow, 1/overflow).


  The bound on x(j) is also used to determine when a step in the


  columnwise method can be performed without fear of overflow.  If


  the computed bound is greater than a large constant, x is scaled to


  prevent overflow, but if the bound overflows, x is set to 0, x(j) to


  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.


  Similarly, a row-wise scheme is used to solve A**T*x = b.  The basic


  algorithm for A upper triangular is


       for j = 1, ..., n


            x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)


       end


  We simultaneously compute two bounds


       G(j) = bound on ( b(i) - A[1:i-1,i]**T * x[1:i-1] ), 1<=i<=j


       M(j) = bound on x(i), 1<=i<=j


  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we


  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.


  Then the bound on x(j) is


       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |


            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )


                      1<=i<=j


  and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater


  than max(underflow, 1/overflow).

 DLATRS3 solves one of the triangular systems


    A * X = B * diag(scale)  or  A**T * X = B * diag(scale)


 with scaling to prevent overflow.  Here A is an upper or lower


 triangular matrix, A**T denotes the transpose of A. X and B are


 n by nrhs matrices and scale is an nrhs element vector of scaling


 factors. A scaling factor scale(j) is usually less than or equal


 to 1, chosen such that X(:,j) is less than the overflow threshold.


 If the matrix A is singular (A(j,j) = 0 for some j), then


 a non-trivial solution to A*X = 0 is returned. If the system is


 so badly scaled that the solution cannot be represented as


 (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned.


 This is a BLAS-3 version of LATRS for solving several right


 hand sides simultaneously.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


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

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

NORMIN

          NORMIN is CHARACTER*1


          Specifies whether CNORM has been set or not.


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


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


                  be computed and stored in CNORM.

N

          N is INTEGER


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

NRHS

          NRHS is INTEGER


          The number of columns of X.  NRHS >= 0.

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          The triangular matrix A.  If UPLO = 'U', the leading n by n


          upper triangular part of the array A contains the upper


          triangular matrix, and the strictly lower triangular part of


          A is not referenced.  If UPLO = 'L', the leading n by n lower


          triangular part of the array A contains the lower triangular


          matrix, and the strictly upper triangular part of A is not


          referenced.  If DIAG = 'U', the diagonal elements of A are


          also not referenced and are assumed to be 1.

LDA

          LDA is INTEGER


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

X

          X is DOUBLE PRECISION array, dimension (LDX,NRHS)


          On entry, the right hand side B of the triangular system.


          On exit, X is overwritten by the solution matrix X.

LDX

          LDX is INTEGER


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

SCALE

          SCALE is DOUBLE PRECISION array, dimension (NRHS)


          The scaling factor s(k) is for the triangular system


          A * x(:,k) = s(k)*b(:,k)  or  A**T* x(:,k) = s(k)*b(:,k).


          If SCALE = 0, the matrix A is singular or badly scaled.


          If A(j,j) = 0 is encountered, a non-trivial vector x(:,k)


          that is an exact or approximate solution to A*x(:,k) = 0


          is returned. If the system so badly scaled that solution


          cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0


          is returned.

CNORM

          CNORM is DOUBLE PRECISION array, dimension (N)


          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)


          contains the norm of the off-diagonal part of the j-th column


          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal


          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)


          must be greater than or equal to the 1-norm.


          If NORMIN = 'N', CNORM is an output argument and CNORM(j)


          returns the 1-norm of the offdiagonal part of the j-th column


          of A.

WORK

          WORK is DOUBLE PRECISION array, dimension (LWORK).


          On exit, if INFO = 0, WORK(1) returns the optimal size of


          WORK.

LWORK LWORK is INTEGER LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAUU2 computes the product U * U**T or L**T * L, where the triangular


 factor U or L is stored in the upper or lower triangular part of


 the array A.


 If UPLO = 'U' or 'u' then the upper triangle of the result is stored,


 overwriting the factor U in A.


 If UPLO = 'L' or 'l' then the lower triangle of the result is stored,


 overwriting the factor L in A.


 This is the unblocked form of the algorithm, calling Level 2 BLAS.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the triangular factor stored in the array A


          is upper or lower triangular:


          = 'U':  Upper triangular


          = 'L':  Lower triangular

N

          N is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the triangular factor U or L.


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


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


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


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

LDA

          LDA is INTEGER


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

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLAUUM computes the product U * U**T or L**T * L, where the triangular


 factor U or L is stored in the upper or lower triangular part of


 the array A.


 If UPLO = 'U' or 'u' then the upper triangle of the result is stored,


 overwriting the factor U in A.


 If UPLO = 'L' or 'l' then the lower triangle of the result is stored,


 overwriting the factor L in A.


 This is the blocked form of the algorithm, calling Level 3 BLAS.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the triangular factor stored in the array A


          is upper or lower triangular:


          = 'U':  Upper triangular


          = 'L':  Lower triangular

N

          N is INTEGER


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

A

          A is DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the triangular factor U or L.


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


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


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


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

LDA

          LDA is INTEGER


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

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DRSCL multiplies an n-element real vector x by the real scalar 1/a.


 This is done without overflow or underflow as long as


 the final result x/a does not overflow or underflow.

N

          N is INTEGER


          The number of components of the vector x.

SA

          SA is DOUBLE PRECISION


          The scalar a which is used to divide each component of x.


          SA must be >= 0, or the subroutine will divide by zero.

SX

          SX is DOUBLE PRECISION array, dimension


                         (1+(N-1)*abs(INCX))


          The n-element vector x.

INCX

          INCX is INTEGER


          The increment between successive values of the vector SX.


          > 0:  SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i),     1< i<= n

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DTPRFB applies a real 'triangular-pentagonal' block reflector H or its


 transpose H**T to a real matrix C, which is composed of two


 blocks A and B, either from the left or right.

SIDE

          SIDE is CHARACTER*1


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


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

TRANS

          TRANS is CHARACTER*1


          = 'N': apply H (No transpose)


          = 'T': apply H**T (Transpose)

DIRECT

          DIRECT is CHARACTER*1


          Indicates how H is formed from a product of elementary


          reflectors


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


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

STOREV

          STOREV is CHARACTER*1


          Indicates how the vectors which define the elementary


          reflectors are stored:


          = 'C': Columns


          = 'R': Rows

M

          M is INTEGER


          The number of rows of the matrix B.


          M >= 0.

N

          N is INTEGER


          The number of columns of the matrix B.


          N >= 0.

K

          K is INTEGER


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


          reflectors whose product defines the block reflector.


          K >= 0.

L

          L is INTEGER


          The order of the trapezoidal part of V.


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

V

          V is DOUBLE PRECISION array, dimension


                                (LDV,K) if STOREV = 'C'


                                (LDV,M) if STOREV = 'R' and SIDE = 'L'


                                (LDV,N) if STOREV = 'R' and SIDE = 'R'


          The pentagonal matrix V, which contains the elementary reflectors


          H(1), H(2), ..., H(K).  See Further Details.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);


          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);


          if STOREV = 'R', LDV >= K.

T

          T is DOUBLE PRECISION array, dimension (LDT,K)


          The triangular K-by-K matrix T in the representation of the


          block reflector.

LDT

          LDT is INTEGER


          The leading dimension of the array T.


          LDT >= K.

A

          A is DOUBLE PRECISION array, dimension


          (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'


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


          On exit, A is overwritten by the corresponding block of


          H*C or H**T*C or C*H or C*H**T.  See Further Details.

LDA

          LDA is INTEGER


          The leading dimension of the array A.


          If SIDE = 'L', LDA >= max(1,K);


          If SIDE = 'R', LDA >= max(1,M).

B

          B is DOUBLE PRECISION array, dimension (LDB,N)


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


          On exit, B is overwritten by the corresponding block of


          H*C or H**T*C or C*H or C*H**T.  See Further Details.

LDB

          LDB is INTEGER


          The leading dimension of the array B.


          LDB >= max(1,M).

WORK

          WORK is DOUBLE PRECISION array, dimension


          (LDWORK,N) if SIDE = 'L',


          (LDWORK,K) if SIDE = 'R'.

LDWORK

          LDWORK is INTEGER


          The leading dimension of the array WORK.


          If SIDE = 'L', LDWORK >= K;


          if SIDE = 'R', LDWORK >= M.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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


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


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


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


                                      [B].


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


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


                                      [A].


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


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


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


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


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


                                         [V2]


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


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


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


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


                                         [V1]


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


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


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


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


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


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


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

 SLATRD reduces NB rows and columns of a real symmetric matrix A to


 symmetric tridiagonal form by an orthogonal similarity


 transformation Q**T * A * Q, and returns the matrices V and W which are


 needed to apply the transformation to the unreduced part of A.


 If UPLO = 'U', SLATRD reduces the last NB rows and columns of a


 matrix, of which the upper triangle is supplied;


 if UPLO = 'L', SLATRD reduces the first NB rows and columns of a


 matrix, of which the lower triangle is supplied.


 This is an auxiliary routine called by SSYTRD.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the upper or lower triangular part of the


          symmetric matrix A is stored:


          = 'U': Upper triangular


          = 'L': Lower triangular

N

          N is INTEGER


          The order of the matrix A.

NB

          NB is INTEGER


          The number of rows and columns to be reduced.

A

          A is REAL array, dimension (LDA,N)


          On entry, the symmetric matrix A.  If UPLO = 'U', the leading


          n-by-n upper triangular part of A contains the upper


          triangular part of the matrix A, and the strictly lower


          triangular part of A is not referenced.  If UPLO = 'L', the


          leading n-by-n lower triangular part of A contains the lower


          triangular part of the matrix A, and the strictly upper


          triangular part of A is not referenced.


          On exit:


          if UPLO = 'U', the last NB columns have been reduced to


            tridiagonal form, with the diagonal elements overwriting


            the diagonal elements of A; the elements above the diagonal


            with the array TAU, represent the orthogonal matrix Q as a


            product of elementary reflectors;


          if UPLO = 'L', the first NB columns have been reduced to


            tridiagonal form, with the diagonal elements overwriting


            the diagonal elements of A; the elements below the diagonal


            with the array TAU, represent the  orthogonal matrix Q as a


            product of elementary reflectors.


          See Further Details.

LDA

          LDA is INTEGER


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

E

          E is REAL array, dimension (N-1)


          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal


          elements of the last NB columns of the reduced matrix;


          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of


          the first NB columns of the reduced matrix.

TAU

          TAU is REAL array, dimension (N-1)


          The scalar factors of the elementary reflectors, stored in


          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.


          See Further Details.

W

          W is REAL array, dimension (LDW,NB)


          The n-by-nb matrix W required to update the unreduced part


          of A.

LDW

          LDW is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


  reflectors


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


  Each H(i) has the form


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


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


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


  and tau in TAU(i-1).


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


  reflectors


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


  Each H(i) has the form


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


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


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


  and tau in TAU(i).


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


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


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


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


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


  with n = 5 and nb = 2:


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


    (  a   a   a   v4  v5 )              (  d                  )


    (      a   a   v4  v5 )              (  1   d              )


    (          a   1   v5 )              (  v1  1   a          )


    (              d   1  )              (  v1  v2  a   a      )


    (                  d  )              (  v1  v2  a   a   a  )


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


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


  an element of the vector defining H(i).

 SLATRS3 solves one of the triangular systems


    A * X = B * diag(scale)  or  A**T * X = B * diag(scale)


 with scaling to prevent overflow.  Here A is an upper or lower


 triangular matrix, A**T denotes the transpose of A. X and B are


 n by nrhs matrices and scale is an nrhs element vector of scaling


 factors. A scaling factor scale(j) is usually less than or equal


 to 1, chosen such that X(:,j) is less than the overflow threshold.


 If the matrix A is singular (A(j,j) = 0 for some j), then


 a non-trivial solution to A*X = 0 is returned. If the system is


 so badly scaled that the solution cannot be represented as


 (1/scale(k))*X(:,k), then x(:,k) = 0 and scale(k) is returned.


 This is a BLAS-3 version of LATRS for solving several right


 hand sides simultaneously.

UPLO

          UPLO is CHARACTER*1


          Specifies whether the matrix A is upper or lower triangular.


          = 'U':  Upper triangular


          = 'L':  Lower triangular

TRANS

          TRANS is CHARACTER*1


          Specifies the operation applied to A.


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


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


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

DIAG

          DIAG is CHARACTER*1


          Specifies whether or not the matrix A is unit triangular.


          = 'N':  Non-unit triangular


          = 'U':  Unit triangular

NORMIN

          NORMIN is CHARACTER*1


          Specifies whether CNORM has been set or not.


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


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


                  be computed and stored in CNORM.

N

          N is INTEGER


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

NRHS

          NRHS is INTEGER


          The number of columns of X.  NRHS >= 0.

A

          A is REAL array, dimension (LDA,N)


          The triangular matrix A.  If UPLO = 'U', the leading n by n


          upper triangular part of the array A contains the upper


          triangular matrix, and the strictly lower triangular part of


          A is not referenced.  If UPLO = 'L', the leading n by n lower


          triangular part of the array A contains the lower triangular


          matrix, and the strictly upper triangular part of A is not


          referenced.  If DIAG = 'U', the diagonal elements of A are


          also not referenced and are assumed to be 1.

LDA

          LDA is INTEGER


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

X

          X is REAL array, dimension (LDX,NRHS)


          On entry, the right hand side B of the triangular system.


          On exit, X is overwritten by the solution matrix X.

LDX

          LDX is INTEGER


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

SCALE

          SCALE is REAL array, dimension (NRHS)


          The scaling factor s(k) is for the triangular system


          A * x(:,k) = s(k)*b(:,k)  or  A**T* x(:,k) = s(k)*b(:,k).


          If SCALE = 0, the matrix A is singular or badly scaled.


          If A(j,j) = 0 is encountered, a non-trivial vector x(:,k)


          that is an exact or approximate solution to A*x(:,k) = 0


          is returned. If the system so badly scaled that solution


          cannot be presented as x(:,k) * 1/s(k), then x(:,k) = 0


          is returned.

CNORM

          CNORM is REAL array, dimension (N)


          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)


          contains the norm of the off-diagonal part of the j-th column


          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal


          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)


          must be greater than or equal to the 1-norm.


          If NORMIN = 'N', CNORM is an output argument and CNORM(j)


          returns the 1-norm of the offdiagonal part of the j-th column


          of A.

WORK

          WORK is REAL array, dimension (LWORK).


          On exit, if INFO = 0, WORK(1) returns the optimal size of


          WORK.

LWORK LWORK is INTEGER LWORK >= MAX(1, 2*NBA * MAX(NBA, MIN(NRHS, 32)), where NBA = (N + NB - 1)/NB and NB is the optimal block size.