SBBCSD computes the CS decomposition of an orthogonal matrix in
 bidiagonal-block form,


     [ B11 | B12 0  0 ]
     [  0  |  0 -I  0 ]
 X = [----------------]
     [ B21 | B22 0  0 ]
     [  0  |  0  0  I ]

                               [  C | -S  0  0 ]
                   [ U1 |    ] [  0 |  0 -I  0 ] [ V1 |    ]**T
                 = [---------] [---------------] [---------]   .
                   [    | U2 ] [  S |  C  0  0 ] [    | V2 ]
                               [  0 |  0  0  I ]

 X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
 than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
 transposed and/or permuted. This can be done in constant time using
 the TRANS and SIGNS options. See SORCSD for details.)

 The bidiagonal matrices B11, B12, B21, and B22 are represented
 implicitly by angles THETA(1:Q) and PHI(1:Q-1).

 The orthogonal matrices U1, U2, V1T, and V2T are input/output.
 The input matrices are pre- or post-multiplied by the appropriate
 singular vector matrices.

          JOBU1 is CHARACTER
          = 'Y':      U1 is updated;
          otherwise:  U1 is not updated.

          JOBU2 is CHARACTER
          = 'Y':      U2 is updated;
          otherwise:  U2 is not updated.

          JOBV1T is CHARACTER
          = 'Y':      V1T is updated;
          otherwise:  V1T is not updated.

          JOBV2T is CHARACTER
          = 'Y':      V2T is updated;
          otherwise:  V2T is not updated.

          TRANS is CHARACTER
          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                      order;
          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                      major order.

          M is INTEGER
          The number of rows and columns in X, the orthogonal matrix in
          bidiagonal-block form.

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

          Q is INTEGER
          The number of columns in the top-left block of X.
          0 <= Q <= MIN(P,M-P,M-Q).

          THETA is REAL array, dimension (Q)
          On entry, the angles THETA(1),...,THETA(Q) that, along with
          PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
          form. On exit, the angles whose cosines and sines define the
          diagonal blocks in the CS decomposition.

          PHI is REAL array, dimension (Q-1)
          The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
          THETA(Q), define the matrix in bidiagonal-block form.

          U1 is REAL array, dimension (LDU1,P)
          On entry, a P-by-P matrix. On exit, U1 is postmultiplied
          by the left singular vector matrix common to [ B11 ; 0 ] and
          [ B12 0 0 ; 0 -I 0 0 ].

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

          U2 is REAL array, dimension (LDU2,M-P)
          On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is
          postmultiplied by the left singular vector matrix common to
          [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].

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

          V1T is REAL array, dimension (LDV1T,Q)
          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
          by the transpose of the right singular vector
          matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].

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

          V2T is REAL array, dimenison (LDV2T,M-Q)
          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
          premultiplied by the transpose of the right
          singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
          [ B22 0 0 ; 0 0 I ].

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

          B11D is REAL array, dimension (Q)
          When SBBCSD converges, B11D contains the cosines of THETA(1),
          ..., THETA(Q). If SBBCSD fails to converge, then B11D
          contains the diagonal of the partially reduced top-left
          block.

          B11E is REAL array, dimension (Q-1)
          When SBBCSD converges, B11E contains zeros. If SBBCSD fails
          to converge, then B11E contains the superdiagonal of the
          partially reduced top-left block.

          B12D is REAL array, dimension (Q)
          When SBBCSD converges, B12D contains the negative sines of
          THETA(1), ..., THETA(Q). If SBBCSD fails to converge, then
          B12D contains the diagonal of the partially reduced top-right
          block.

          B12E is REAL array, dimension (Q-1)
          When SBBCSD converges, B12E contains zeros. If SBBCSD fails
          to converge, then B12E contains the subdiagonal of the
          partially reduced top-right block.

          B21D is REAL array, dimension (Q)
          When SBBCSD converges, B21D contains the negative sines of
          THETA(1), ..., THETA(Q). If SBBCSD fails to converge, then
          B21D contains the diagonal of the partially reduced bottom-left
          block.

          B21E is REAL array, dimension (Q-1)
          When SBBCSD converges, B21E contains zeros. If SBBCSD fails
          to converge, then B21E contains the subdiagonal of the
          partially reduced bottom-left block.

          B22D is REAL array, dimension (Q)
          When SBBCSD converges, B22D contains the negative sines of
          THETA(1), ..., THETA(Q). If SBBCSD fails to converge, then
          B22D contains the diagonal of the partially reduced bottom-right
          block.

          B22E is REAL array, dimension (Q-1)
          When SBBCSD converges, B22E contains zeros. If SBBCSD fails
          to converge, then B22E contains the subdiagonal of the
          partially reduced bottom-right block.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= MAX(1,8*Q).

          If LWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal size of the WORK array,
          returns this value as the first entry of the work array, and
          no error message related to LWORK is issued by XERBLA.

          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if SBBCSD did not converge, INFO specifies the number
                of nonzero entries in PHI, and B11D, B11E, etc.,
                contain the partially reduced matrix.

  TOLMUL  REAL, default = MAX(10,MIN(100,EPS**(-1/8)))
          TOLMUL controls the convergence criterion of the QR loop.
          Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
          are within TOLMUL*EPS of either bound.

 SGGHD3 reduces a pair of real matrices (A,B) to generalized upper
 Hessenberg form using orthogonal transformations, where A is a
 general matrix and B is upper triangular.  The form of the
 generalized eigenvalue problem is
    A*x = lambda*B*x,
 and B is typically made upper triangular by computing its QR
 factorization and moving the orthogonal matrix Q to the left side
 of the equation.

 This subroutine simultaneously reduces A to a Hessenberg matrix H:
    Q**T*A*Z = H
 and transforms B to another upper triangular matrix T:
    Q**T*B*Z = T
 in order to reduce the problem to its standard form
    H*y = lambda*T*y
 where y = Z**T*x.

 The orthogonal matrices Q and Z are determined as products of Givens
 rotations.  They may either be formed explicitly, or they may be
 postmultiplied into input matrices Q1 and Z1, so that

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

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

 If Q1 is the orthogonal matrix from the QR factorization of B in the
 original equation A*x = lambda*B*x, then SGGHD3 reduces the original
 problem to generalized Hessenberg form.

 This is a blocked variant of SGGHRD, using matrix-matrix
 multiplications for parts of the computation to enhance performance.

          COMPQ is CHARACTER*1
          = 'N': do not compute Q;
          = 'I': Q is initialized to the unit matrix, and the
                 orthogonal matrix Q is returned;
          = 'V': Q must contain an orthogonal matrix Q1 on entry,
                 and the product Q1*Q is returned.

          COMPZ is CHARACTER*1
          = 'N': do not compute Z;
          = 'I': Z is initialized to the unit matrix, and the
                 orthogonal matrix Z is returned;
          = 'V': Z must contain an orthogonal matrix Z1 on entry,
                 and the product Z1*Z is returned.

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

          ILO is INTEGER

          IHI is INTEGER

          ILO and IHI mark the rows and columns of A which are to be
          reduced.  It is assumed that A is already upper triangular
          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
          normally set by a previous call to SGGBAL; otherwise they
          should be set to 1 and N respectively.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

          A is REAL array, dimension (LDA, N)
          On entry, the N-by-N general matrix to be reduced.
          On exit, the upper triangle and the first subdiagonal of A
          are overwritten with the upper Hessenberg matrix H, and the
          rest is set to zero.

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

          B is REAL array, dimension (LDB, N)
          On entry, the N-by-N upper triangular matrix B.
          On exit, the upper triangular matrix T = Q**T B Z.  The
          elements below the diagonal are set to zero.

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

          Q is REAL array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the orthogonal matrix Q1,
          typically from the QR factorization of B.
          On exit, if COMPQ='I', the orthogonal matrix Q, and if
          COMPQ = 'V', the product Q1*Q.
          Not referenced if COMPQ='N'.

          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.

          Z is REAL array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the orthogonal matrix Z1.
          On exit, if COMPZ='I', the orthogonal matrix Z, and if
          COMPZ = 'V', the product Z1*Z.
          Not referenced if COMPZ='N'.

          LDZ is INTEGER
          The leading dimension of the array Z.
          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.

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

          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 1.
          For optimum performance LWORK >= 6*N*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

  This routine reduces A to Hessenberg form and maintains B in
  using a blocked variant of Moler and Stewart's original algorithm,
  as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
  (BIT 2008).

 SGGHRD reduces a pair of real matrices (A,B) to generalized upper
 Hessenberg form using orthogonal transformations, where A is a
 general matrix and B is upper triangular.  The form of the
 generalized eigenvalue problem is
    A*x = lambda*B*x,
 and B is typically made upper triangular by computing its QR
 factorization and moving the orthogonal matrix Q to the left side
 of the equation.

 This subroutine simultaneously reduces A to a Hessenberg matrix H:
    Q**T*A*Z = H
 and transforms B to another upper triangular matrix T:
    Q**T*B*Z = T
 in order to reduce the problem to its standard form
    H*y = lambda*T*y
 where y = Z**T*x.

 The orthogonal matrices Q and Z are determined as products of Givens
 rotations.  They may either be formed explicitly, or they may be
 postmultiplied into input matrices Q1 and Z1, so that

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

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

 If Q1 is the orthogonal matrix from the QR factorization of B in the
 original equation A*x = lambda*B*x, then SGGHRD reduces the original
 problem to generalized Hessenberg form.

          COMPQ is CHARACTER*1
          = 'N': do not compute Q;
          = 'I': Q is initialized to the unit matrix, and the
                 orthogonal matrix Q is returned;
          = 'V': Q must contain an orthogonal matrix Q1 on entry,
                 and the product Q1*Q is returned.

          COMPZ is CHARACTER*1
          = 'N': do not compute Z;
          = 'I': Z is initialized to the unit matrix, and the
                 orthogonal matrix Z is returned;
          = 'V': Z must contain an orthogonal matrix Z1 on entry,
                 and the product Z1*Z is returned.

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

          ILO is INTEGER

          IHI is INTEGER

          ILO and IHI mark the rows and columns of A which are to be
          reduced.  It is assumed that A is already upper triangular
          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
          normally set by a previous call to SGGBAL; otherwise they
          should be set to 1 and N respectively.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

          A is REAL array, dimension (LDA, N)
          On entry, the N-by-N general matrix to be reduced.
          On exit, the upper triangle and the first subdiagonal of A
          are overwritten with the upper Hessenberg matrix H, and the
          rest is set to zero.

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

          B is REAL array, dimension (LDB, N)
          On entry, the N-by-N upper triangular matrix B.
          On exit, the upper triangular matrix T = Q**T B Z.  The
          elements below the diagonal are set to zero.

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

          Q is REAL array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the orthogonal matrix Q1,
          typically from the QR factorization of B.
          On exit, if COMPQ='I', the orthogonal matrix Q, and if
          COMPQ = 'V', the product Q1*Q.
          Not referenced if COMPQ='N'.

          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.

          Z is REAL array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the orthogonal matrix Z1.
          On exit, if COMPZ='I', the orthogonal matrix Z, and if
          COMPZ = 'V', the product Z1*Z.
          Not referenced if COMPZ='N'.

          LDZ is INTEGER
          The leading dimension of the array Z.
          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.

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

  This routine reduces A to Hessenberg and B to triangular form by
  an unblocked reduction, as described in _Matrix_Computations_,
  by Golub and Van Loan (Johns Hopkins Press.)

 SGGQRF computes a generalized QR factorization of an N-by-M matrix A
 and an N-by-P matrix B:

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

 where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
 matrix, and R and T assume one of the forms:

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

 where R11 is upper triangular, and

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

 where T12 or T21 is upper triangular.

 In particular, if B is square and nonsingular, the GQR factorization
 of A and B implicitly gives the QR factorization of inv(B)*A:

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

 where inv(B) denotes the inverse of the matrix B, and Z**T denotes the
 transpose of the matrix Z.

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

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

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

          A is REAL array, dimension (LDA,M)
          On entry, the N-by-M matrix A.
          On exit, the elements on and above the diagonal of the array
          contain the min(N,M)-by-M upper trapezoidal matrix R (R is
          upper triangular if N >= M); the elements below the diagonal,
          with the array TAUA, represent the orthogonal matrix Q as a
          product of min(N,M) elementary reflectors (see Further
          Details).

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

          TAUA is REAL array, dimension (min(N,M))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix Q (see Further Details).

          B is REAL array, dimension (LDB,P)
          On entry, the N-by-P matrix B.
          On exit, if N <= P, the upper triangle of the subarray
          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
          if N > P, the elements on and above the (N-P)-th subdiagonal
          contain the N-by-P upper trapezoidal matrix T; the remaining
          elements, with the array TAUB, represent the orthogonal
          matrix Z as a product of elementary reflectors (see Further
          Details).

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

          TAUB is REAL array, dimension (min(N,P))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix Z (see Further Details).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,N,M,P).
          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
          where NB1 is the optimal blocksize for the QR factorization
          of an N-by-M matrix, NB2 is the optimal blocksize for the
          RQ factorization of an N-by-P matrix, and NB3 is the optimal
          blocksize for a call of SORMQR.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

  The matrix Q is represented as a product of elementary reflectors

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

  Each H(i) has the form

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

  where taua is a real scalar, and v is a real vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
  and taua in TAUA(i).
  To form Q explicitly, use LAPACK subroutine SORGQR.
  To use Q to update another matrix, use LAPACK subroutine SORMQR.

  The matrix Z is represented as a product of elementary reflectors

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

  Each H(i) has the form

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

  where taub is a real scalar, and v is a real vector with
  v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
  B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
  To form Z explicitly, use LAPACK subroutine SORGRQ.
  To use Z to update another matrix, use LAPACK subroutine SORMRQ.

 SGGRQF computes a generalized RQ factorization of an M-by-N matrix A
 and a P-by-N matrix B:

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

 where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
 matrix, and R and T assume one of the forms:

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

 where R12 or R21 is upper triangular, and

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

 where T11 is upper triangular.

 In particular, if B is square and nonsingular, the GRQ factorization
 of A and B implicitly gives the RQ factorization of A*inv(B):

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

 where inv(B) denotes the inverse of the matrix B, and Z**T denotes the
 transpose of the matrix Z.

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

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

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

          A is REAL array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, if M <= N, the upper triangle of the subarray
          A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R;
          if M > N, the elements on and above the (M-N)-th subdiagonal
          contain the M-by-N upper trapezoidal matrix R; the remaining
          elements, with the array TAUA, represent the orthogonal
          matrix Q as a product of elementary reflectors (see Further
          Details).

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

          TAUA is REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix Q (see Further Details).

          B is REAL array, dimension (LDB,N)
          On entry, the P-by-N matrix B.
          On exit, the elements on and above the diagonal of the array
          contain the min(P,N)-by-N upper trapezoidal matrix T (T is
          upper triangular if P >= N); the elements below the diagonal,
          with the array TAUB, represent the orthogonal matrix Z as a
          product of elementary reflectors (see Further Details).

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

          TAUB is REAL array, dimension (min(P,N))
          The scalar factors of the elementary reflectors which
          represent the orthogonal matrix Z (see Further Details).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,N,M,P).
          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
          where NB1 is the optimal blocksize for the RQ factorization
          of an M-by-N matrix, NB2 is the optimal blocksize for the
          QR factorization of a P-by-N matrix, and NB3 is the optimal
          blocksize for a call of SORMRQ.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INF0= -i, the i-th argument had an illegal value.

  The matrix Q is represented as a product of elementary reflectors

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

  Each H(i) has the form

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

  where taua is a real scalar, and v is a real vector with
  v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
  A(m-k+i,1:n-k+i-1), and taua in TAUA(i).
  To form Q explicitly, use LAPACK subroutine SORGRQ.
  To use Q to update another matrix, use LAPACK subroutine SORMRQ.

  The matrix Z is represented as a product of elementary reflectors

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

  Each H(i) has the form

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

  where taub is a real scalar, and v is a real vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i),
  and taub in TAUB(i).
  To form Z explicitly, use LAPACK subroutine SORGQR.
  To use Z to update another matrix, use LAPACK subroutine SORMQR.

 SGGSVP3 computes orthogonal matrices U, V and Q such that

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

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

                  N-K-L  K    L
  V**T*B*Q =   L ( 0     0   B13 )
             P-L ( 0     0    0  )

 where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
 upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
 otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
 numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T.

 This decomposition is the preprocessing step for computing the
 Generalized Singular Value Decomposition (GSVD), see subroutine
 SGGSVD3.

          JOBU is CHARACTER*1
          = 'U':  Orthogonal matrix U is computed;
          = 'N':  U is not computed.

          JOBV is CHARACTER*1
          = 'V':  Orthogonal matrix V is computed;
          = 'N':  V is not computed.

          JOBQ is CHARACTER*1
          = 'Q':  Orthogonal matrix Q is computed;
          = 'N':  Q is not computed.

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

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

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

          A is REAL array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, A contains the triangular (or trapezoidal) matrix
          described in the Purpose section.

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

          B is REAL array, dimension (LDB,N)
          On entry, the P-by-N matrix B.
          On exit, B contains the triangular matrix described in
          the Purpose section.

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

          TOLA is REAL

          TOLB is REAL

          TOLA and TOLB are the thresholds to determine the effective
          numerical rank of matrix B and a subblock of A. Generally,
          they are set to
             TOLA = MAX(M,N)*norm(A)*MACHEPS,
             TOLB = MAX(P,N)*norm(B)*MACHEPS.
          The size of TOLA and TOLB may affect the size of backward
          errors of the decomposition.

          K is INTEGER

          L is INTEGER

          On exit, K and L specify the dimension of the subblocks
          described in Purpose section.
          K + L = effective numerical rank of (A**T,B**T)**T.

          U is REAL array, dimension (LDU,M)
          If JOBU = 'U', U contains the orthogonal matrix U.
          If JOBU = 'N', U is not referenced.

          LDU is INTEGER
          The leading dimension of the array U. LDU >= max(1,M) if
          JOBU = 'U'; LDU >= 1 otherwise.

          V is REAL array, dimension (LDV,P)
          If JOBV = 'V', V contains the orthogonal matrix V.
          If JOBV = 'N', V is not referenced.

          LDV is INTEGER
          The leading dimension of the array V. LDV >= max(1,P) if
          JOBV = 'V'; LDV >= 1 otherwise.

          Q is REAL array, dimension (LDQ,N)
          If JOBQ = 'Q', Q contains the orthogonal matrix Q.
          If JOBQ = 'N', Q is not referenced.

          LDQ is INTEGER
          The leading dimension of the array Q. LDQ >= max(1,N) if
          JOBQ = 'Q'; LDQ >= 1 otherwise.

          IWORK is INTEGER array, dimension (N)

          TAU is REAL array, dimension (N)

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

  The subroutine uses LAPACK subroutine SGEQP3 for the QR factorization
  with column pivoting to detect the effective numerical rank of the
  a matrix. It may be replaced by a better rank determination strategy.

  SGGSVP3 replaces the deprecated subroutine SGGSVP.

 SGSVJ0 is called from SGESVJ as a pre-processor and that is its main
 purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
 it does not check convergence (stopping criterion). Few tuning
 parameters (marked by [TP]) are available for the implementer.

          JOBV is CHARACTER*1
          Specifies whether the output from this procedure is used
          to compute the matrix V:
          = 'V': the product of the Jacobi rotations is accumulated
                 by postmulyiplying the N-by-N array V.
                (See the description of V.)
          = 'A': the product of the Jacobi rotations is accumulated
                 by postmulyiplying the MV-by-N array V.
                (See the descriptions of MV and V.)
          = 'N': the Jacobi rotations are not accumulated.

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

          N is INTEGER
          The number of columns of the input matrix A.
          M >= N >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, M-by-N matrix A, such that A*diag(D) represents
          the input matrix.
          On exit,
          A_onexit * D_onexit represents the input matrix A*diag(D)
          post-multiplied by a sequence of Jacobi rotations, where the
          rotation threshold and the total number of sweeps are given in
          TOL and NSWEEP, respectively.
          (See the descriptions of D, TOL and NSWEEP.)

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

          D is REAL array, dimension (N)
          The array D accumulates the scaling factors from the fast scaled
          Jacobi rotations.
          On entry, A*diag(D) represents the input matrix.
          On exit, A_onexit*diag(D_onexit) represents the input matrix
          post-multiplied by a sequence of Jacobi rotations, where the
          rotation threshold and the total number of sweeps are given in
          TOL and NSWEEP, respectively.
          (See the descriptions of A, TOL and NSWEEP.)

          SVA is REAL array, dimension (N)
          On entry, SVA contains the Euclidean norms of the columns of
          the matrix A*diag(D).
          On exit, SVA contains the Euclidean norms of the columns of
          the matrix onexit*diag(D_onexit).

          MV is INTEGER
          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV = 'N',   then MV is not referenced.

          V is REAL array, dimension (LDV,N)
          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV = 'N',   then V is not referenced.

          LDV is INTEGER
          The leading dimension of the array V,  LDV >= 1.
          If JOBV = 'V', LDV .GE. N.
          If JOBV = 'A', LDV .GE. MV.

          EPS is REAL
          EPS = SLAMCH('Epsilon')

          SFMIN is REAL
          SFMIN = SLAMCH('Safe Minimum')

          TOL is REAL
          TOL is the threshold for Jacobi rotations. For a pair
          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.

          NSWEEP is INTEGER
          NSWEEP is the number of sweeps of Jacobi rotations to be
          performed.

          WORK is REAL array, dimension LWORK.

          LWORK is INTEGER
          LWORK is the dimension of WORK. LWORK .GE. M.

          INFO is INTEGER
          = 0 : successful exit.
          < 0 : if INFO = -i, then the i-th argument had an illegal value

 SGSVJ1 is called from SGESVJ as a pre-processor and that is its main
 purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
 it targets only particular pivots and it does not check convergence
 (stopping criterion). Few tunning parameters (marked by [TP]) are
 available for the implementer.

 Further Details
 ~~~~~~~~~~~~~~~
 SGSVJ1 applies few sweeps of Jacobi rotations in the column space of
 the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
 off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
 block-entries (tiles) of the (1,2) off-diagonal block are marked by the
 [x]'s in the following scheme:

    | *  *  * [x] [x] [x]|
    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
    |[x] [x] [x] *  *  * |
    |[x] [x] [x] *  *  * |
    |[x] [x] [x] *  *  * |

 In terms of the columns of A, the first N1 columns are rotated 'against'
 the remaining N-N1 columns, trying to increase the angle between the
 corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
 tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
 The number of sweeps is given in NSWEEP and the orthogonality threshold
 is given in TOL.

          JOBV is CHARACTER*1
          Specifies whether the output from this procedure is used
          to compute the matrix V:
          = 'V': the product of the Jacobi rotations is accumulated
                 by postmulyiplying the N-by-N array V.
                (See the description of V.)
          = 'A': the product of the Jacobi rotations is accumulated
                 by postmulyiplying the MV-by-N array V.
                (See the descriptions of MV and V.)
          = 'N': the Jacobi rotations are not accumulated.

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

          N is INTEGER
          The number of columns of the input matrix A.
          M >= N >= 0.

          N1 is INTEGER
          N1 specifies the 2 x 2 block partition, the first N1 columns are
          rotated 'against' the remaining N-N1 columns of A.

          A is REAL array, dimension (LDA,N)
          On entry, M-by-N matrix A, such that A*diag(D) represents
          the input matrix.
          On exit,
          A_onexit * D_onexit represents the input matrix A*diag(D)
          post-multiplied by a sequence of Jacobi rotations, where the
          rotation threshold and the total number of sweeps are given in
          TOL and NSWEEP, respectively.
          (See the descriptions of N1, D, TOL and NSWEEP.)

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

          D is REAL array, dimension (N)
          The array D accumulates the scaling factors from the fast scaled
          Jacobi rotations.
          On entry, A*diag(D) represents the input matrix.
          On exit, A_onexit*diag(D_onexit) represents the input matrix
          post-multiplied by a sequence of Jacobi rotations, where the
          rotation threshold and the total number of sweeps are given in
          TOL and NSWEEP, respectively.
          (See the descriptions of N1, A, TOL and NSWEEP.)

          SVA is REAL array, dimension (N)
          On entry, SVA contains the Euclidean norms of the columns of
          the matrix A*diag(D).
          On exit, SVA contains the Euclidean norms of the columns of
          the matrix onexit*diag(D_onexit).

          MV is INTEGER
          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV = 'N',   then MV is not referenced.

          V is REAL array, dimension (LDV,N)
          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV = 'N',   then V is not referenced.

          LDV is INTEGER
          The leading dimension of the array V,  LDV >= 1.
          If JOBV = 'V', LDV .GE. N.
          If JOBV = 'A', LDV .GE. MV.

          EPS is REAL
          EPS = SLAMCH('Epsilon')

          SFMIN is REAL
          SFMIN = SLAMCH('Safe Minimum')

          TOL is REAL
          TOL is the threshold for Jacobi rotations. For a pair
          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.

          NSWEEP is INTEGER
          NSWEEP is the number of sweeps of Jacobi rotations to be
          performed.

          WORK is REAL array, dimension LWORK.

          LWORK is INTEGER
          LWORK is the dimension of WORK. LWORK .GE. M.

          INFO is INTEGER
          = 0 : successful exit.
          < 0 : if INFO = -i, then the i-th argument had an illegal value

 SHSEIN uses inverse iteration to find specified right and/or left
 eigenvectors of a real upper Hessenberg matrix H.

 The right eigenvector x and the left eigenvector y of the matrix H
 corresponding to an eigenvalue w are defined by:

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

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

          SIDE is CHARACTER*1
          = 'R': compute right eigenvectors only;
          = 'L': compute left eigenvectors only;
          = 'B': compute both right and left eigenvectors.

          EIGSRC is CHARACTER*1
          Specifies the source of eigenvalues supplied in (WR,WI):
          = 'Q': the eigenvalues were found using SHSEQR; thus, if
                 H has zero subdiagonal elements, and so is
                 block-triangular, then the j-th eigenvalue can be
                 assumed to be an eigenvalue of the block containing
                 the j-th row/column.  This property allows SHSEIN to
                 perform inverse iteration on just one diagonal block.
          = 'N': no assumptions are made on the correspondence
                 between eigenvalues and diagonal blocks.  In this
                 case, SHSEIN must always perform inverse iteration
                 using the whole matrix H.

          INITV is CHARACTER*1
          = 'N': no initial vectors are supplied;
          = 'U': user-supplied initial vectors are stored in the arrays
                 VL and/or VR.

          SELECT is LOGICAL array, dimension (N)
          Specifies the eigenvectors to be computed. To select the
          real eigenvector corresponding to a real eigenvalue WR(j),
          SELECT(j) must be set to .TRUE.. To select the complex
          eigenvector corresponding to a complex eigenvalue
          (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
          either SELECT(j) or SELECT(j+1) or both must be set to
          .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
          .FALSE..

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

          H is REAL array, dimension (LDH,N)
          The upper Hessenberg matrix H.
          If a NaN is detected in H, the routine will return with INFO=-6.

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

          WR is REAL array, dimension (N)

          WI is REAL array, dimension (N)

          On entry, the real and imaginary parts of the eigenvalues of
          H; a complex conjugate pair of eigenvalues must be stored in
          consecutive elements of WR and WI.
          On exit, WR may have been altered since close eigenvalues
          are perturbed slightly in searching for independent
          eigenvectors.

          VL is REAL array, dimension (LDVL,MM)
          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
          contain starting vectors for the inverse iteration for the
          left eigenvectors; the starting vector for each eigenvector
          must be in the same column(s) in which the eigenvector will
          be stored.
          On exit, if SIDE = 'L' or 'B', the left eigenvectors
          specified by SELECT will be stored consecutively in the
          columns of VL, in the same order as their eigenvalues. A
          complex eigenvector corresponding to a complex eigenvalue is
          stored in two consecutive columns, the first holding the real
          part and the second the imaginary part.
          If SIDE = 'R', VL is not referenced.

          LDVL is INTEGER
          The leading dimension of the array VL.
          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.

          VR is REAL array, dimension (LDVR,MM)
          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
          contain starting vectors for the inverse iteration for the
          right eigenvectors; the starting vector for each eigenvector
          must be in the same column(s) in which the eigenvector will
          be stored.
          On exit, if SIDE = 'R' or 'B', the right eigenvectors
          specified by SELECT will be stored consecutively in the
          columns of VR, in the same order as their eigenvalues. A
          complex eigenvector corresponding to a complex eigenvalue is
          stored in two consecutive columns, the first holding the real
          part and the second the imaginary part.
          If SIDE = 'L', VR is not referenced.

          LDVR is INTEGER
          The leading dimension of the array VR.
          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.

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

          M is INTEGER
          The number of columns in the arrays VL and/or VR required to
          store the eigenvectors; each selected real eigenvector
          occupies one column and each selected complex eigenvector
          occupies two columns.

          WORK is REAL array, dimension ((N+2)*N)

          IFAILL is INTEGER array, dimension (MM)
          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
          eigenvector in the i-th column of VL (corresponding to the
          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
          eigenvector converged satisfactorily. If the i-th and (i+1)th
          columns of VL hold a complex eigenvector, then IFAILL(i) and
          IFAILL(i+1) are set to the same value.
          If SIDE = 'R', IFAILL is not referenced.

          IFAILR is INTEGER array, dimension (MM)
          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
          eigenvector in the i-th column of VR (corresponding to the
          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
          eigenvector converged satisfactorily. If the i-th and (i+1)th
          columns of VR hold a complex eigenvector, then IFAILR(i) and
          IFAILR(i+1) are set to the same value.
          If SIDE = 'L', IFAILR is not referenced.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, i is the number of eigenvectors which
                failed to converge; see IFAILL and IFAILR for further
                details.

  Each eigenvector is normalized so that the element of largest
  magnitude has magnitude 1; here the magnitude of a complex number
  (x,y) is taken to be |x|+|y|.

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

          JOB is CHARACTER*1
           = 'E':  compute eigenvalues only;
           = 'S':  compute eigenvalues and the Schur form T.

          COMPZ is CHARACTER*1
           = 'N':  no Schur vectors are computed;
           = 'I':  Z is initialized to the unit matrix and the matrix Z
                   of Schur vectors of H is returned;
           = 'V':  Z must contain an orthogonal matrix Q on entry, and
                   the product Q*Z is returned.

          N is INTEGER
           The order of the matrix H.  N .GE. 0.

          ILO is INTEGER

          IHI is INTEGER

           It is assumed that H is already upper triangular in rows
           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
           set by a previous call to SGEBAL, and then passed to ZGEHRD
           when the matrix output by SGEBAL is reduced to Hessenberg
           form. Otherwise ILO and IHI should be set to 1 and N
           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
           If N = 0, then ILO = 1 and IHI = 0.

          H is REAL array, dimension (LDH,N)
           On entry, the upper Hessenberg matrix H.
           On exit, if INFO = 0 and JOB = 'S', 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).LT.0. If INFO = 0 and JOB = 'E', the
           contents of H are unspecified on exit.  (The output value of
           H when INFO.GT.0 is given under the description of INFO
           below.)

           Unlike earlier versions of SHSEQR, this subroutine may
           explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
           or j = IHI+1, IHI+2, ... N.

          LDH is INTEGER
           The leading dimension of the array H. LDH .GE. max(1,N).

          WR is REAL array, dimension (N)

          WI is REAL array, dimension (N)

           The real and imaginary parts, respectively, of the computed
           eigenvalues. 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) .GT. 0 and
           WI(i+1) .LT. 0. If JOB = 'S', 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).

          Z is REAL array, dimension (LDZ,N)
           If COMPZ = 'N', Z is not referenced.
           If COMPZ = 'I', on entry Z need not be set and on exit,
           if INFO = 0, Z contains the orthogonal matrix Z of the Schur
           vectors of H.  If COMPZ = 'V', on entry Z must contain an
           N-by-N matrix Q, which is assumed to be equal to the unit
           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
           if INFO = 0, Z contains Q*Z.
           Normally Q is the orthogonal matrix generated by SORGHR
           after the call to SGEHRD which formed the Hessenberg matrix
           H. (The output value of Z when INFO.GT.0 is given under
           the description of INFO below.)

          LDZ is INTEGER
           The leading dimension of the array Z.  if COMPZ = 'I' or
           COMPZ = 'V', then LDZ.GE.MAX(1,N).  Otherwize, LDZ.GE.1.

          WORK is REAL array, dimension (LWORK)
           On exit, if INFO = 0, WORK(1) returns an estimate of
           the optimal value for LWORK.

          LWORK is INTEGER
           The dimension of the array WORK.  LWORK .GE. max(1,N)
           is sufficient and delivers very good and sometimes
           optimal performance.  However, LWORK as large as 11*N
           may be required for optimal performance.  A workspace
           query is recommended to determine the optimal workspace
           size.

           If LWORK = -1, then SHSEQR does a workspace query.
           In this case, SHSEQR 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 is INTEGER
             =  0:  successful exit
           .LT. 0:  if INFO = -i, the i-th argument had an illegal
                    value
           .GT. 0:  if INFO = i, SHSEQR failed to compute all of
                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
                and WI contain those eigenvalues which have been
                successfully computed.  (Failures are rare.)

                If INFO .GT. 0 and JOB = 'E', then on exit, the
                remaining unconverged eigenvalues are the eigen-
                values of the upper Hessenberg matrix rows and
                columns ILO through INFO of the final, output
                value of H.

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

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

                where U is the orthogonal matrix in (*) (regard-
                less of the value of JOB.)

                If INFO .GT. 0 and COMPZ = 'I', then on exit
                      (final value of Z)  = U
                where U is the orthogonal matrix in (*) (regard-
                less of the value of JOB.)

                If INFO .GT. 0 and COMPZ = 'N', then Z is not
                accessed.

             Default values supplied by
             ILAENV(ISPEC,'SHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
             It is suggested that these defaults be adjusted in order
             to attain best performance in each particular
             computational environment.

            ISPEC=12: The SLAHQR vs SLAQR0 crossover point.
                      Default: 75. (Must be at least 11.)

            ISPEC=13: Recommended deflation window size.
                      This depends on ILO, IHI and NS.  NS is the
                      number of simultaneous shifts returned
                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)
                      The default for (IHI-ILO+1).LE.500 is NS.
                      The default for (IHI-ILO+1).GT.500 is 3*NS/2.

            ISPEC=14: Nibble crossover point. (See IPARMQ for
                      details.)  Default: 14% of deflation window
                      size.

            ISPEC=15: Number of simultaneous shifts in a multishift
                      QR iteration.

                      If IHI-ILO+1 is ...

                      greater than      ...but less    ... the
                      or equal to ...      than        default is

                           1               30          NS =   2(+)
                          30               60          NS =   4(+)
                          60              150          NS =  10(+)
                         150              590          NS =  **
                         590             3000          NS =  64
                        3000             6000          NS = 128
                        6000             infinity      NS = 256

                  (+)  By default some or all matrices of this order
                       are passed to the implicit double shift routine
                       SLAHQR and this parameter is ignored.  See
                       ISPEC=12 above and comments in IPARMQ for
                       details.

                 (**)  The asterisks (**) indicate an ad-hoc
                       function of N increasing from 10 to 64.

            ISPEC=16: Select structured matrix multiply.
                      If the number of simultaneous shifts (specified
                      by ISPEC=15) is less than 14, then the default
                      for ISPEC=16 is 0.  Otherwise the default for
                      ISPEC=16 is 2.

    SLA_LIN_BERR computes componentwise relative backward error from
    the formula
        max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
    where abs(Z) is the componentwise absolute value of the matrix
    or vector Z.

          N is INTEGER
     The number of linear equations, i.e., the order of the
     matrix A.  N >= 0.

          NZ is INTEGER
     We add (NZ+1)*SLAMCH( 'Safe minimum' ) to R(i) in the numerator to
     guard against spuriously zero residuals. Default value is N.

          NRHS is INTEGER
     The number of right hand sides, i.e., the number of columns
     of the matrices AYB, RES, and BERR.  NRHS >= 0.

          RES is REAL array, dimension (N,NRHS)
     The residual matrix, i.e., the matrix R in the relative backward
     error formula above.

          AYB is REAL array, dimension (N, NRHS)
     The denominator in the relative backward error formula above, i.e.,
     the matrix abs(op(A_s))*abs(Y) + abs(B_s). The matrices A, Y, and B
     are from iterative refinement (see sla_gerfsx_extended.f).

          BERR is REAL array, dimension (NRHS)
     The componentwise relative backward error from the formula above.

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

    This works for all extant IBM's hex and binary floating point
    arithmetics, but not for decimal.

          N is INTEGER
            The length of vectors X, Y, and W.

          X is REAL array, dimension (N)
            The first part of the doubled-single accumulation vector.

          Y is REAL array, dimension (N)
            The second part of the doubled-single accumulation vector.

          W is REAL array, dimension (N)
            The vector to be added.

 SLALS0 applies back the multiplying factors of either the left or the
 right singular vector matrix of a diagonal matrix appended by a row
 to the right hand side matrix B in solving the least squares problem
 using the divide-and-conquer SVD approach.

 For the left singular vector matrix, three types of orthogonal
 matrices are involved:

 (1L) Givens rotations: the number of such rotations is GIVPTR; the
      pairs of columns/rows they were applied to are stored in GIVCOL;
      and the C- and S-values of these rotations are stored in GIVNUM.

 (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
      row, and for J=2:N, PERM(J)-th row of B is to be moved to the
      J-th row.

 (3L) The left singular vector matrix of the remaining matrix.

 For the right singular vector matrix, four types of orthogonal
 matrices are involved:

 (1R) The right singular vector matrix of the remaining matrix.

 (2R) If SQRE = 1, one extra Givens rotation to generate the right
      null space.

 (3R) The inverse transformation of (2L).

 (4R) The inverse transformation of (1L).

          ICOMPQ is INTEGER
         Specifies whether singular vectors are to be computed in
         factored form:
         = 0: Left singular vector matrix.
         = 1: Right singular vector matrix.

          NL is INTEGER
         The row dimension of the upper block. NL >= 1.

          NR is INTEGER
         The row dimension of the lower block. NR >= 1.

          SQRE is INTEGER
         = 0: the lower block is an NR-by-NR square matrix.
         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

         The bidiagonal matrix has row dimension N = NL + NR + 1,
         and column dimension M = N + SQRE.

          NRHS is INTEGER
         The number of columns of B and BX. NRHS must be at least 1.

          B is REAL array, dimension ( LDB, NRHS )
         On input, B contains the right hand sides of the least
         squares problem in rows 1 through M. On output, B contains
         the solution X in rows 1 through N.

          LDB is INTEGER
         The leading dimension of B. LDB must be at least
         max(1,MAX( M, N ) ).

          BX is REAL array, dimension ( LDBX, NRHS )

          LDBX is INTEGER
         The leading dimension of BX.

          PERM is INTEGER array, dimension ( N )
         The permutations (from deflation and sorting) applied
         to the two blocks.

          GIVPTR is INTEGER
         The number of Givens rotations which took place in this
         subproblem.

          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
         Each pair of numbers indicates a pair of rows/columns
         involved in a Givens rotation.

          LDGCOL is INTEGER
         The leading dimension of GIVCOL, must be at least N.

          GIVNUM is REAL array, dimension ( LDGNUM, 2 )
         Each number indicates the C or S value used in the
         corresponding Givens rotation.

          LDGNUM is INTEGER
         The leading dimension of arrays DIFR, POLES and
         GIVNUM, must be at least K.

          POLES is REAL array, dimension ( LDGNUM, 2 )
         On entry, POLES(1:K, 1) contains the new singular
         values obtained from solving the secular equation, and
         POLES(1:K, 2) is an array containing the poles in the secular
         equation.

          DIFL is REAL array, dimension ( K ).
         On entry, DIFL(I) is the distance between I-th updated
         (undeflated) singular value and the I-th (undeflated) old
         singular value.

          DIFR is REAL array, dimension ( LDGNUM, 2 ).
         On entry, DIFR(I, 1) contains the distances between I-th
         updated (undeflated) singular value and the I+1-th
         (undeflated) old singular value. And DIFR(I, 2) is the
         normalizing factor for the I-th right singular vector.

          Z is REAL array, dimension ( K )
         Contain the components of the deflation-adjusted updating row
         vector.

          K is INTEGER
         Contains the dimension of the non-deflated matrix,
         This is the order of the related secular equation. 1 <= K <=N.

          C is REAL
         C contains garbage if SQRE =0 and the C-value of a Givens
         rotation related to the right null space if SQRE = 1.

          S is REAL
         S contains garbage if SQRE =0 and the S-value of a Givens
         rotation related to the right null space if SQRE = 1.

          WORK is REAL array, dimension ( K )

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

 SLALSA is an itermediate step in solving the least squares problem
 by computing the SVD of the coefficient matrix in compact form (The
 singular vectors are computed as products of simple orthorgonal
 matrices.).

 If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector
 matrix of an upper bidiagonal matrix to the right hand side; and if
 ICOMPQ = 1, SLALSA applies the right singular vector matrix to the
 right hand side. The singular vector matrices were generated in
 compact form by SLALSA.

          ICOMPQ is INTEGER
         Specifies whether the left or the right singular vector
         matrix is involved.
         = 0: Left singular vector matrix
         = 1: Right singular vector matrix

          SMLSIZ is INTEGER
         The maximum size of the subproblems at the bottom of the
         computation tree.

          N is INTEGER
         The row and column dimensions of the upper bidiagonal matrix.

          NRHS is INTEGER
         The number of columns of B and BX. NRHS must be at least 1.

          B is REAL array, dimension ( LDB, NRHS )
         On input, B contains the right hand sides of the least
         squares problem in rows 1 through M.
         On output, B contains the solution X in rows 1 through N.

          LDB is INTEGER
         The leading dimension of B in the calling subprogram.
         LDB must be at least max(1,MAX( M, N ) ).

          BX is REAL array, dimension ( LDBX, NRHS )
         On exit, the result of applying the left or right singular
         vector matrix to B.

          LDBX is INTEGER
         The leading dimension of BX.

          U is REAL array, dimension ( LDU, SMLSIZ ).
         On entry, U contains the left singular vector matrices of all
         subproblems at the bottom level.

          LDU is INTEGER, LDU = > N.
         The leading dimension of arrays U, VT, DIFL, DIFR,
         POLES, GIVNUM, and Z.

          VT is REAL array, dimension ( LDU, SMLSIZ+1 ).
         On entry, VT**T contains the right singular vector matrices of
         all subproblems at the bottom level.

          K is INTEGER array, dimension ( N ).

          DIFL is REAL array, dimension ( LDU, NLVL ).
         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.

          DIFR is REAL array, dimension ( LDU, 2 * NLVL ).
         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
         distances between singular values on the I-th level and
         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
         record the normalizing factors of the right singular vectors
         matrices of subproblems on I-th level.

          Z is REAL array, dimension ( LDU, NLVL ).
         On entry, Z(1, I) contains the components of the deflation-
         adjusted updating row vector for subproblems on the I-th
         level.

          POLES is REAL array, dimension ( LDU, 2 * NLVL ).
         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
         singular values involved in the secular equations on the I-th
         level.

          GIVPTR is INTEGER array, dimension ( N ).
         On entry, GIVPTR( I ) records the number of Givens
         rotations performed on the I-th problem on the computation
         tree.

          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
         locations of Givens rotations performed on the I-th level on
         the computation tree.

          LDGCOL is INTEGER, LDGCOL = > N.
         The leading dimension of arrays GIVCOL and PERM.

          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
         On entry, PERM(*, I) records permutations done on the I-th
         level of the computation tree.

          GIVNUM is REAL array, dimension ( LDU, 2 * NLVL ).
         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
         values of Givens rotations performed on the I-th level on the
         computation tree.

          C is REAL array, dimension ( N ).
         On entry, if the I-th subproblem is not square,
         C( I ) contains the C-value of a Givens rotation related to
         the right null space of the I-th subproblem.

          S is REAL array, dimension ( N ).
         On entry, if the I-th subproblem is not square,
         S( I ) contains the S-value of a Givens rotation related to
         the right null space of the I-th subproblem.

          WORK is REAL array.
         The dimension must be at least N.

          IWORK is INTEGER array.
         The dimension must be at least 3 * N

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

 SLALSD uses the singular value decomposition of A to solve the least
 squares problem of finding X to minimize the Euclidean norm of each
 column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
 are N-by-NRHS. The solution X overwrites B.

 The singular values of A smaller than RCOND times the largest
 singular value are treated as zero in solving the least squares
 problem; in this case a minimum norm solution is returned.
 The actual singular values are returned in D in ascending order.

 This code makes very mild assumptions about floating point
 arithmetic. It will work on machines with a guard digit in
 add/subtract, or on those binary machines without guard digits
 which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.

          UPLO is CHARACTER*1
         = 'U': D and E define an upper bidiagonal matrix.
         = 'L': D and E define a  lower bidiagonal matrix.

          SMLSIZ is INTEGER
         The maximum size of the subproblems at the bottom of the
         computation tree.

          N is INTEGER
         The dimension of the  bidiagonal matrix.  N >= 0.

          NRHS is INTEGER
         The number of columns of B. NRHS must be at least 1.

          D is REAL array, dimension (N)
         On entry D contains the main diagonal of the bidiagonal
         matrix. On exit, if INFO = 0, D contains its singular values.

          E is REAL array, dimension (N-1)
         Contains the super-diagonal entries of the bidiagonal matrix.
         On exit, E has been destroyed.

          B is REAL array, dimension (LDB,NRHS)
         On input, B contains the right hand sides of the least
         squares problem. On output, B contains the solution X.

          LDB is INTEGER
         The leading dimension of B in the calling subprogram.
         LDB must be at least max(1,N).

          RCOND is REAL
         The singular values of A less than or equal to RCOND times
         the largest singular value are treated as zero in solving
         the least squares problem. If RCOND is negative,
         machine precision is used instead.
         For example, if diag(S)*X=B were the least squares problem,
         where diag(S) is a diagonal matrix of singular values, the
         solution would be X(i) = B(i) / S(i) if S(i) is greater than
         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
         RCOND*max(S).

          RANK is INTEGER
         The number of singular values of A greater than RCOND times
         the largest singular value.

          WORK is REAL array, dimension at least
         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).

          IWORK is INTEGER array, dimension at least
         (3*N*NLVL + 11*N)

          INFO is INTEGER
         = 0:  successful exit.
         < 0:  if INFO = -i, the i-th argument had an illegal value.
         > 0:  The algorithm failed to compute a singular value while
               working on the submatrix lying in rows and columns
               INFO/(N+1) through MOD(INFO,N+1).

 SLANSF 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 in RFP format.

    SLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
             (
             ( norm1(A),         NORM = '1', 'O' or 'o'
             (
             ( normI(A),         NORM = 'I' or 'i'
             (
             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

 where  norm1  denotes the  one norm of a matrix (maximum column sum),
 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
 normF  denotes the  Frobenius norm of a matrix (square root of sum of
 squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.

          NORM is CHARACTER*1
          Specifies the value to be returned in SLANSF as described
          above.

          TRANSR is CHARACTER*1
          Specifies whether the RFP format of A is normal or
          transposed format.
          = 'N':  RFP format is Normal;
          = 'T':  RFP format is Transpose.

          UPLO is CHARACTER*1
           On entry, UPLO specifies whether the RFP matrix A came from
           an upper or lower triangular matrix as follows:
           = 'U': RFP A came from an upper triangular matrix;
           = 'L': RFP A came from a lower triangular matrix.

          N is INTEGER
          The order of the matrix A. N >= 0. When N = 0, SLANSF is
          set to zero.

          A is REAL array, dimension ( N*(N+1)/2 );
          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
          part of the symmetric matrix A stored in RFP format. See the
          "Notes" below for more details.
          Unchanged on exit.

          WORK is REAL array, dimension (MAX(1,LWORK)),
          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
          WORK is not referenced.

  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

 SLARSCL2 performs a reciprocal diagonal scaling on an vector:
   x <-- inv(D) * x
 where the diagonal matrix D is stored as a vector.

 Eventually to be replaced by BLAS_sge_diag_scale in the new BLAS
 standard.

          M is INTEGER
     The number of rows of D and X. M >= 0.

          N is INTEGER
     The number of columns of X. N >= 0.

          D is REAL array, length M
     Diagonal matrix D, stored as a vector of length M.

          X is REAL array, dimension (LDX,N)
     On entry, the vector X to be scaled by D.
     On exit, the scaled vector.

          LDX is INTEGER
     The leading dimension of the vector X. LDX >= M.

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


 H is a product of k elementary reflectors as returned by STZRZF.

          SIDE is CHARACTER*1
          = 'L': form  H * C
          = 'R': form  C * H

          M is INTEGER
          The number of rows of the matrix C.

          N is INTEGER
          The number of columns of the matrix C.

          L is INTEGER
          The number of entries of the vector V containing
          the meaningful part of the Householder vectors.
          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

          V is REAL array, dimension (1+(L-1)*abs(INCV))
          The vector v in the representation of H as returned by
          STZRZF. V is not used if TAU = 0.

          INCV is INTEGER
          The increment between elements of v. INCV <> 0.

          TAU is REAL
          The value tau in the representation of H.

          C is REAL 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 is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                         (N) if SIDE = 'L'
                      or (M) if SIDE = 'R'

 SLARZB applies a real block reflector H or its transpose H**T to
 a real distributed M-by-N  C from the left or the right.

 Currently, only STOREV = 'R' and DIRECT = 'B' are supported.

          SIDE is CHARACTER*1
          = 'L': apply H or H**T from the Left
          = 'R': apply H or H**T from the Right

          TRANS is CHARACTER*1
          = 'N': apply H (No transpose)
          = 'C': apply H**T (Transpose)

          DIRECT is CHARACTER*1
          Indicates how H is formed from a product of elementary
          reflectors
          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
          = 'B': H = H(k) . . . H(2) H(1) (Backward)

          STOREV is CHARACTER*1
          Indicates how the vectors which define the elementary
          reflectors are stored:
          = 'C': Columnwise                        (not supported yet)
          = 'R': Rowwise

          M is INTEGER
          The number of rows of the matrix C.

          N is INTEGER
          The number of columns of the matrix C.

          K is INTEGER
          The order of the matrix T (= the number of elementary
          reflectors whose product defines the block reflector).

          L is INTEGER
          The number of columns of the matrix V containing the
          meaningful part of the Householder reflectors.
          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

          V is REAL array, dimension (LDV,NV).
          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.

          LDV is INTEGER
          The leading dimension of the array V.
          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.

          T is REAL array, dimension (LDT,K)
          The triangular K-by-K matrix T in the representation of the
          block reflector.

          LDT is INTEGER
          The leading dimension of the array T. LDT >= K.

          C is REAL 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 is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (LDWORK,K)

          LDWORK is INTEGER
          The leading dimension of the array WORK.
          If SIDE = 'L', LDWORK >= max(1,N);
          if SIDE = 'R', LDWORK >= max(1,M).

 SLARZT 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

 Currently, only STOREV = 'R' and DIRECT = 'B' are supported.

          DIRECT is CHARACTER*1
          Specifies the order in which the elementary reflectors are
          multiplied to form the block reflector:
          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
          = 'B': H = H(k) . . . H(2) H(1) (Backward)

          STOREV is CHARACTER*1
          Specifies how the vectors which define the elementary
          reflectors are stored (see also Further Details):
          = 'C': columnwise                        (not supported yet)
          = 'R': rowwise

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

          K is INTEGER
          The order of the triangular factor T (= the number of
          elementary reflectors). K >= 1.

          V is REAL array, dimension
                               (LDV,K) if STOREV = 'C'
                               (LDV,N) if STOREV = 'R'
          The matrix V. See further details.

          LDV is INTEGER
          The leading dimension of the array V.
          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i).

          T is REAL 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 is INTEGER
          The leading dimension of the array T. LDT >= K.

  The shape of the matrix V and the storage of the vectors which define
  the H(i) is best illustrated by the following example with n = 5 and
  k = 3. The elements equal to 1 are not stored; the corresponding
  array elements are modified but restored on exit. The rest of the
  array is not used.

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

                                              ______V_____
         ( v1 v2 v3 )                        /                     ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )
     V = ( v1 v2 v3 )                      ( v2 v2 v2 v2 v2 . . . 1   )
         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )
         ( v1 v2 v3 )
            .  .  .
            .  .  .
            1  .  .
               1  .
                  1

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

                                                        ______V_____
            1                                          /                        .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )
            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )
            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )
            .  .  .
         ( v1 v2 v3 )
         ( v1 v2 v3 )
     V = ( v1 v2 v3 )
         ( v1 v2 v3 )
         ( v1 v2 v3 )

 SLASCL2 performs a diagonal scaling on a vector:
   x <-- D * x
 where the diagonal matrix D is stored as a vector.

 Eventually to be replaced by BLAS_sge_diag_scale in the new BLAS
 standard.

          M is INTEGER
     The number of rows of D and X. M >= 0.

          N is INTEGER
     The number of columns of X. N >= 0.

          D is REAL array, length M
     Diagonal matrix D, stored as a vector of length M.

          X is REAL array, dimension (LDX,N)
     On entry, the vector X to be scaled by D.
     On exit, the scaled vector.

          LDX is INTEGER
     The leading dimension of the vector X. LDX >= M.

 SLATRZ factors the M-by-(M+L) real upper trapezoidal matrix
 [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R  0 ) * Z, by means
 of orthogonal transformations.  Z is an (M+L)-by-(M+L) orthogonal
 matrix and, R and A1 are M-by-M upper triangular matrices.

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

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

          L is INTEGER
          The number of columns of the matrix A containing the
          meaningful part of the Householder vectors. N-M >= L >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the leading M-by-N upper trapezoidal part of the
          array A must contain the matrix to be factorized.
          On exit, the leading M-by-M upper triangular part of A
          contains the upper triangular matrix R, and elements N-L+1 to
          N of the first M rows of A, with the array TAU, represent the
          orthogonal matrix Z as a product of M elementary reflectors.

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

          TAU is REAL array, dimension (M)
          The scalar factors of the elementary reflectors.

          WORK is REAL array, dimension (M)

  The factorization is obtained by Householder's method.  The kth
  transformation matrix, Z( k ), which is used to introduce zeros into
  the ( m - k + 1 )th row of A, is given in the form

     Z( k ) = ( I     0   ),
              ( 0  T( k ) )

  where

     T( k ) = I - tau*u( k )*u( k )**T,   u( k ) = (   1    ),
                                                 (   0    )
                                                 ( z( k ) )

  tau is a scalar and z( k ) is an l element vector. tau and z( k )
  are chosen to annihilate the elements of the kth row of A2.

  The scalar tau is returned in the kth element of TAU and the vector
  u( k ) in the kth row of A2, such that the elements of z( k ) are
  in  a( k, l + 1 ), ..., a( k, n ). The elements of R are returned in
  the upper triangular part of A1.

  Z is given by

     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).

 SOPGTR generates a real orthogonal matrix Q which is defined as the
 product of n-1 elementary reflectors H(i) of order n, as returned by
 SSPTRD using packed storage:

 if UPLO = 'U', Q = H(n-1) . . . H(2) H(1),

 if UPLO = 'L', Q = H(1) H(2) . . . H(n-1).

          UPLO is CHARACTER*1
          = 'U': Upper triangular packed storage used in previous
                 call to SSPTRD;
          = 'L': Lower triangular packed storage used in previous
                 call to SSPTRD.

          N is INTEGER
          The order of the matrix Q. N >= 0.

          AP is REAL array, dimension (N*(N+1)/2)
          The vectors which define the elementary reflectors, as
          returned by SSPTRD.

          TAU is REAL array, dimension (N-1)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SSPTRD.

          Q is REAL array, dimension (LDQ,N)
          The N-by-N orthogonal matrix Q.

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

          WORK is REAL array, dimension (N-1)

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

 SOPMTR overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix of order nq, with nq = m if
 SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
 nq-1 elementary reflectors, as returned by SSPTRD using packed
 storage:

 if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);

 if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          UPLO is CHARACTER*1
          = 'U': Upper triangular packed storage used in previous
                 call to SSPTRD;
          = 'L': Lower triangular packed storage used in previous
                 call to SSPTRD.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          AP is REAL array, dimension
                               (M*(M+1)/2) if SIDE = 'L'
                               (N*(N+1)/2) if SIDE = 'R'
          The vectors which define the elementary reflectors, as
          returned by SSPTRD.  AP is modified by the routine but
          restored on exit.

          TAU is REAL array, dimension (M-1) if SIDE = 'L'
                                     or (N-1) if SIDE = 'R'
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SSPTRD.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                                   (N) if SIDE = 'L'
                                   (M) if SIDE = 'R'

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

 SORBDB simultaneously bidiagonalizes the blocks of an M-by-M
 partitioned orthogonal matrix X:

                                 [ B11 | B12 0  0 ]
     [ X11 | X12 ]   [ P1 |    ] [  0  |  0 -I  0 ] [ Q1 |    ]**T
 X = [-----------] = [---------] [----------------] [---------]   .
     [ X21 | X22 ]   [    | P2 ] [ B21 | B22 0  0 ] [    | Q2 ]
                                 [  0  |  0  0  I ]

 X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
 not the case, then X must be transposed and/or permuted. This can be
 done in constant time using the TRANS and SIGNS options. See SORCSD
 for details.)

 The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
 (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
 represented implicitly by Householder vectors.

 B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
 implicitly by angles THETA, PHI.

          TRANS is CHARACTER
          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                      order;
          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                      major order.

          SIGNS is CHARACTER
          = 'O':      The lower-left block is made nonpositive (the
                      "other" convention);
          otherwise:  The upper-right block is made nonpositive (the
                      "default" convention).

          M is INTEGER
          The number of rows and columns in X.

          P is INTEGER
          The number of rows in X11 and X12. 0 <= P <= M.

          Q is INTEGER
          The number of columns in X11 and X21. 0 <= Q <=
          MIN(P,M-P,M-Q).

          X11 is REAL array, dimension (LDX11,Q)
          On entry, the top-left block of the orthogonal matrix to be
          reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the columns of tril(X11) specify reflectors for P1,
             the rows of triu(X11,1) specify reflectors for Q1;
          else TRANS = 'T', and
             the rows of triu(X11) specify reflectors for P1,
             the columns of tril(X11,-1) specify reflectors for Q1.

          LDX11 is INTEGER
          The leading dimension of X11. If TRANS = 'N', then LDX11 >=
          P; else LDX11 >= Q.

          X12 is REAL array, dimension (LDX12,M-Q)
          On entry, the top-right block of the orthogonal matrix to
          be reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the rows of triu(X12) specify the first P reflectors for
             Q2;
          else TRANS = 'T', and
             the columns of tril(X12) specify the first P reflectors
             for Q2.

          LDX12 is INTEGER
          The leading dimension of X12. If TRANS = 'N', then LDX12 >=
          P; else LDX11 >= M-Q.

          X21 is REAL array, dimension (LDX21,Q)
          On entry, the bottom-left block of the orthogonal matrix to
          be reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the columns of tril(X21) specify reflectors for P2;
          else TRANS = 'T', and
             the rows of triu(X21) specify reflectors for P2.

          LDX21 is INTEGER
          The leading dimension of X21. If TRANS = 'N', then LDX21 >=
          M-P; else LDX21 >= Q.

          X22 is REAL array, dimension (LDX22,M-Q)
          On entry, the bottom-right block of the orthogonal matrix to
          be reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
             M-P-Q reflectors for Q2,
          else TRANS = 'T', and
             the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
             M-P-Q reflectors for P2.

          LDX22 is INTEGER
          The leading dimension of X22. If TRANS = 'N', then LDX22 >=
          M-P; else LDX22 >= M-Q.

          THETA is REAL array, dimension (Q)
          The entries of the bidiagonal blocks B11, B12, B21, B22 can
          be computed from the angles THETA and PHI. See Further
          Details.

          PHI is REAL array, dimension (Q-1)
          The entries of the bidiagonal blocks B11, B12, B21, B22 can
          be computed from the angles THETA and PHI. See Further
          Details.

          TAUP1 is REAL array, dimension (P)
          The scalar factors of the elementary reflectors that define
          P1.

          TAUP2 is REAL array, dimension (M-P)
          The scalar factors of the elementary reflectors that define
          P2.

          TAUQ1 is REAL array, dimension (Q)
          The scalar factors of the elementary reflectors that define
          Q1.

          TAUQ2 is REAL array, dimension (M-Q)
          The scalar factors of the elementary reflectors that define
          Q2.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= M-Q.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

  The bidiagonal blocks B11, B12, B21, and B22 are represented
  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
  lower bidiagonal. Every entry in each bidiagonal band is a product
  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
  [1] or SORCSD for details.

  P1, P2, Q1, and Q2 are represented as products of elementary
  reflectors. See SORCSD for details on generating P1, P2, Q1, and Q2
  using SORGQR and SORGLQ.

M 

          M is INTEGER
           The number of rows X11 plus the number of rows in X21.

          P is INTEGER
           The number of rows in X11. 0 <= P <= M.

          Q is INTEGER
           The number of columns in X11 and X21. 0 <= Q <=
           MIN(P,M-P,M-Q).

          X11 is REAL array, dimension (LDX11,Q)
           On entry, the top block of the matrix X to be reduced. On
           exit, the columns of tril(X11) specify reflectors for P1 and
           the rows of triu(X11,1) specify reflectors for Q1.

          LDX11 is INTEGER
           The leading dimension of X11. LDX11 >= P.

          X21 is REAL array, dimension (LDX21,Q)
           On entry, the bottom block of the matrix X to be reduced. On
           exit, the columns of tril(X21) specify reflectors for P2.

          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= M-P.

          THETA is REAL array, dimension (Q)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          PHI is REAL array, dimension (Q-1)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          TAUP1 is REAL array, dimension (P)
           The scalar factors of the elementary reflectors that define
           P1.

          TAUP2 is REAL array, dimension (M-P)
           The scalar factors of the elementary reflectors that define
           P2.

          TAUQ1 is REAL array, dimension (Q)
           The scalar factors of the elementary reflectors that define
           Q1.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= M-Q.

           If LWORK = -1, then a workspace query is assumed; the routine
           only calculates the optimal size of the WORK array, returns
           this value as the first entry of the WORK array, and no error
           message related to LWORK is issued by XERBLA.

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

  The upper-bidiagonal blocks B11, B21 are represented implicitly by
  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
  in each bidiagonal band is a product of a sine or cosine of a THETA
  with a sine or cosine of a PHI. See [1] or SORCSD for details.

  P1, P2, and Q1 are represented as products of elementary reflectors.
  See SORCSD2BY1 for details on generating P1, P2, and Q1 using SORGQR
  and SORGLQ.

M 

          M is INTEGER
           The number of rows X11 plus the number of rows in X21.

          P is INTEGER
           The number of rows in X11. 0 <= P <= min(M-P,Q,M-Q).

          Q is INTEGER
           The number of columns in X11 and X21. 0 <= Q <= M.

          X11 is REAL array, dimension (LDX11,Q)
           On entry, the top block of the matrix X to be reduced. On
           exit, the columns of tril(X11) specify reflectors for P1 and
           the rows of triu(X11,1) specify reflectors for Q1.

          LDX11 is INTEGER
           The leading dimension of X11. LDX11 >= P.

          X21 is REAL array, dimension (LDX21,Q)
           On entry, the bottom block of the matrix X to be reduced. On
           exit, the columns of tril(X21) specify reflectors for P2.

          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= M-P.

          THETA is REAL array, dimension (Q)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          PHI is REAL array, dimension (Q-1)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          TAUP1 is REAL array, dimension (P)
           The scalar factors of the elementary reflectors that define
           P1.

          TAUP2 is REAL array, dimension (M-P)
           The scalar factors of the elementary reflectors that define
           P2.

          TAUQ1 is REAL array, dimension (Q)
           The scalar factors of the elementary reflectors that define
           Q1.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= M-Q.

           If LWORK = -1, then a workspace query is assumed; the routine
           only calculates the optimal size of the WORK array, returns
           this value as the first entry of the WORK array, and no error
           message related to LWORK is issued by XERBLA.

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

  The upper-bidiagonal blocks B11, B21 are represented implicitly by
  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
  in each bidiagonal band is a product of a sine or cosine of a THETA
  with a sine or cosine of a PHI. See [1] or SORCSD for details.

  P1, P2, and Q1 are represented as products of elementary reflectors.
  See SORCSD2BY1 for details on generating P1, P2, and Q1 using SORGQR
  and SORGLQ.

M 

          M is INTEGER
           The number of rows X11 plus the number of rows in X21.

          P is INTEGER
           The number of rows in X11. 0 <= P <= M. M-P <= min(P,Q,M-Q).

          Q is INTEGER
           The number of columns in X11 and X21. 0 <= Q <= M.

          X11 is REAL array, dimension (LDX11,Q)
           On entry, the top block of the matrix X to be reduced. On
           exit, the columns of tril(X11) specify reflectors for P1 and
           the rows of triu(X11,1) specify reflectors for Q1.

          LDX11 is INTEGER
           The leading dimension of X11. LDX11 >= P.

          X21 is REAL array, dimension (LDX21,Q)
           On entry, the bottom block of the matrix X to be reduced. On
           exit, the columns of tril(X21) specify reflectors for P2.

          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= M-P.

          THETA is REAL array, dimension (Q)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          PHI is REAL array, dimension (Q-1)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          TAUP1 is REAL array, dimension (P)
           The scalar factors of the elementary reflectors that define
           P1.

          TAUP2 is REAL array, dimension (M-P)
           The scalar factors of the elementary reflectors that define
           P2.

          TAUQ1 is REAL array, dimension (Q)
           The scalar factors of the elementary reflectors that define
           Q1.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= M-Q.

           If LWORK = -1, then a workspace query is assumed; the routine
           only calculates the optimal size of the WORK array, returns
           this value as the first entry of the WORK array, and no error
           message related to LWORK is issued by XERBLA.

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

  The upper-bidiagonal blocks B11, B21 are represented implicitly by
  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
  in each bidiagonal band is a product of a sine or cosine of a THETA
  with a sine or cosine of a PHI. See [1] or SORCSD for details.

  P1, P2, and Q1 are represented as products of elementary reflectors.
  See SORCSD2BY1 for details on generating P1, P2, and Q1 using SORGQR
  and SORGLQ.

M 

          M is INTEGER
           The number of rows X11 plus the number of rows in X21.

          P is INTEGER
           The number of rows in X11. 0 <= P <= M.

          Q is INTEGER
           The number of columns in X11 and X21. 0 <= Q <= M and
           M-Q <= min(P,M-P,Q).

          X11 is REAL array, dimension (LDX11,Q)
           On entry, the top block of the matrix X to be reduced. On
           exit, the columns of tril(X11) specify reflectors for P1 and
           the rows of triu(X11,1) specify reflectors for Q1.

          LDX11 is INTEGER
           The leading dimension of X11. LDX11 >= P.

          X21 is REAL array, dimension (LDX21,Q)
           On entry, the bottom block of the matrix X to be reduced. On
           exit, the columns of tril(X21) specify reflectors for P2.

          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= M-P.

          THETA is REAL array, dimension (Q)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          PHI is REAL array, dimension (Q-1)
           The entries of the bidiagonal blocks B11, B21 are defined by
           THETA and PHI. See Further Details.

          TAUP1 is REAL array, dimension (P)
           The scalar factors of the elementary reflectors that define
           P1.

          TAUP2 is REAL array, dimension (M-P)
           The scalar factors of the elementary reflectors that define
           P2.

          TAUQ1 is REAL array, dimension (Q)
           The scalar factors of the elementary reflectors that define
           Q1.

          PHANTOM is REAL array, dimension (M)
           The routine computes an M-by-1 column vector Y that is
           orthogonal to the columns of [ X11; X21 ]. PHANTOM(1:P) and
           PHANTOM(P+1:M) contain Householder vectors for Y(1:P) and
           Y(P+1:M), respectively.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= M-Q.

           If LWORK = -1, then a workspace query is assumed; the routine
           only calculates the optimal size of the WORK array, returns
           this value as the first entry of the WORK array, and no error
           message related to LWORK is issued by XERBLA.

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

  The upper-bidiagonal blocks B11, B21 are represented implicitly by
  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
  in each bidiagonal band is a product of a sine or cosine of a THETA
  with a sine or cosine of a PHI. See [1] or SORCSD for details.

  P1, P2, and Q1 are represented as products of elementary reflectors.
  See SORCSD2BY1 for details on generating P1, P2, and Q1 using SORGQR
  and SORGLQ.

M1 

          M1 is INTEGER
           The dimension of X1 and the number of rows in Q1. 0 <= M1.

          M2 is INTEGER
           The dimension of X2 and the number of rows in Q2. 0 <= M2.

          N is INTEGER
           The number of columns in Q1 and Q2. 0 <= N.

          X1 is REAL array, dimension (M1)
           On entry, the top part of the vector to be orthogonalized.
           On exit, the top part of the projected vector.

          INCX1 is INTEGER
           Increment for entries of X1.

          X2 is REAL array, dimension (M2)
           On entry, the bottom part of the vector to be
           orthogonalized. On exit, the bottom part of the projected
           vector.

          INCX2 is INTEGER
           Increment for entries of X2.

          Q1 is REAL array, dimension (LDQ1, N)
           The top part of the orthonormal basis matrix.

          LDQ1 is INTEGER
           The leading dimension of Q1. LDQ1 >= M1.

          Q2 is REAL array, dimension (LDQ2, N)
           The bottom part of the orthonormal basis matrix.

          LDQ2 is INTEGER
           The leading dimension of Q2. LDQ2 >= M2.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= N.

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

M1 

          M1 is INTEGER
           The dimension of X1 and the number of rows in Q1. 0 <= M1.

          M2 is INTEGER
           The dimension of X2 and the number of rows in Q2. 0 <= M2.

          N is INTEGER
           The number of columns in Q1 and Q2. 0 <= N.

          X1 is REAL array, dimension (M1)
           On entry, the top part of the vector to be orthogonalized.
           On exit, the top part of the projected vector.

          INCX1 is INTEGER
           Increment for entries of X1.

          X2 is REAL array, dimension (M2)
           On entry, the bottom part of the vector to be
           orthogonalized. On exit, the bottom part of the projected
           vector.

          INCX2 is INTEGER
           Increment for entries of X2.

          Q1 is REAL array, dimension (LDQ1, N)
           The top part of the orthonormal basis matrix.

          LDQ1 is INTEGER
           The leading dimension of Q1. LDQ1 >= M1.

          Q2 is REAL array, dimension (LDQ2, N)
           The bottom part of the orthonormal basis matrix.

          LDQ2 is INTEGER
           The leading dimension of Q2. LDQ2 >= M2.

          WORK is REAL array, dimension (LWORK)

          LWORK is INTEGER
           The dimension of the array WORK. LWORK >= N.

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

 SORCSD computes the CS decomposition of an M-by-M partitioned
 orthogonal matrix X:

                                 [  I  0  0 |  0  0  0 ]
                                 [  0  C  0 |  0 -S  0 ]
     [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**T
 X = [-----------] = [---------] [---------------------] [---------]   .
     [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ]
                                 [  0  S  0 |  0  C  0 ]
                                 [  0  0  I |  0  0  0 ]

 X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
 (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
 R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
 which R = MIN(P,M-P,Q,M-Q).

          JOBU1 is CHARACTER
          = 'Y':      U1 is computed;
          otherwise:  U1 is not computed.

          JOBU2 is CHARACTER
          = 'Y':      U2 is computed;
          otherwise:  U2 is not computed.

          JOBV1T is CHARACTER
          = 'Y':      V1T is computed;
          otherwise:  V1T is not computed.

          JOBV2T is CHARACTER
          = 'Y':      V2T is computed;
          otherwise:  V2T is not computed.

          TRANS is CHARACTER
          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                      order;
          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                      major order.

          SIGNS is CHARACTER
          = 'O':      The lower-left block is made nonpositive (the
                      "other" convention);
          otherwise:  The upper-right block is made nonpositive (the
                      "default" convention).

          M is INTEGER
          The number of rows and columns in X.

          P is INTEGER
          The number of rows in X11 and X12. 0 <= P <= M.

          Q is INTEGER
          The number of columns in X11 and X21. 0 <= Q <= M.

          X11 is REAL array, dimension (LDX11,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.

          LDX11 is INTEGER
          The leading dimension of X11. LDX11 >= MAX(1,P).

          X12 is REAL array, dimension (LDX12,M-Q)
          On entry, part of the orthogonal matrix whose CSD is desired.

          LDX12 is INTEGER
          The leading dimension of X12. LDX12 >= MAX(1,P).

          X21 is REAL array, dimension (LDX21,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.

          LDX21 is INTEGER
          The leading dimension of X11. LDX21 >= MAX(1,M-P).

          X22 is REAL array, dimension (LDX22,M-Q)
          On entry, part of the orthogonal matrix whose CSD is desired.

          LDX22 is INTEGER
          The leading dimension of X11. LDX22 >= MAX(1,M-P).

          THETA is REAL array, dimension (R), in which R =
          MIN(P,M-P,Q,M-Q).
          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).

          U1 is REAL array, dimension (P)
          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.

          LDU1 is INTEGER
          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
          MAX(1,P).

          U2 is REAL array, dimension (M-P)
          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
          matrix U2.

          LDU2 is INTEGER
          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
          MAX(1,M-P).

          V1T is REAL array, dimension (Q)
          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
          matrix V1**T.

          LDV1T is INTEGER
          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
          MAX(1,Q).

          V2T is REAL array, dimension (M-Q)
          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal
          matrix V2**T.

          LDV2T is INTEGER
          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
          MAX(1,M-Q).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
          define the matrix in intermediate bidiagonal-block form
          remaining after nonconvergence. INFO specifies the number
          of nonzero PHI's.

          LWORK is INTEGER
          The dimension of the array WORK.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the work array, and no error
          message related to LWORK is issued by XERBLA.

          IWORK is INTEGER array, dimension (M-MIN(P, M-P, Q, M-Q))

          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  SBBCSD did not converge. See the description of WORK
                above for details.

JOBU1 

          JOBU1 is CHARACTER
          = 'Y':      U1 is computed;
          otherwise:  U1 is not computed.

          JOBU2 is CHARACTER
          = 'Y':      U2 is computed;
          otherwise:  U2 is not computed.

          JOBV1T is CHARACTER
          = 'Y':      V1T is computed;
          otherwise:  V1T is not computed.

          M is INTEGER
          The number of rows in X.

          P is INTEGER
          The number of rows in X11. 0 <= P <= M.

          Q is INTEGER
          The number of columns in X11 and X21. 0 <= Q <= M.

          X11 is REAL array, dimension (LDX11,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.

          LDX11 is INTEGER
          The leading dimension of X11. LDX11 >= MAX(1,P).

          X21 is REAL array, dimension (LDX21,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.

          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= MAX(1,M-P).

          THETA is REAL array, dimension (R), in which R =
          MIN(P,M-P,Q,M-Q).
          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).

          U1 is REAL array, dimension (P)
          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.

          LDU1 is INTEGER
          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
          MAX(1,P).

          U2 is REAL array, dimension (M-P)
          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
          matrix U2.

          LDU2 is INTEGER
          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
          MAX(1,M-P).

          V1T is REAL array, dimension (Q)
          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
          matrix V1**T.

          LDV1T is INTEGER
          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
          MAX(1,Q).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
          define the matrix in intermediate bidiagonal-block form
          remaining after nonconvergence. INFO specifies the number
          of nonzero PHI's.

          LWORK is INTEGER
          The dimension of the array WORK.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the work array, and no error
          message related to LWORK is issued by XERBLA.

          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))

          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  SBBCSD did not converge. See the description of WORK
                above for details.

 SORG2L generates an m by n real matrix Q with orthonormal columns,
 which is defined as the last n columns of a product of k elementary
 reflectors of order m

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

 as returned by SGEQLF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. M >= N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. N >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the (n-k+i)-th column must contain the vector which
          defines the elementary reflector H(i), for i = 1,2,...,k, as
          returned by SGEQLF in the last k columns of its array
          argument A.
          On exit, the m by n matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQLF.

          WORK is REAL array, dimension (N)

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument has an illegal value

 SORG2R generates an m by n real matrix Q with orthonormal columns,
 which is defined as the first n columns of a product of k elementary
 reflectors of order m

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

 as returned by SGEQRF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. M >= N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. N >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the i-th column must contain the vector which
          defines the elementary reflector H(i), for i = 1,2,...,k, as
          returned by SGEQRF in the first k columns of its array
          argument A.
          On exit, the m-by-n matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQRF.

          WORK is REAL array, dimension (N)

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument has an illegal value

 SORGHR generates a real orthogonal matrix Q which is defined as the
 product of IHI-ILO elementary reflectors of order N, as returned by
 SGEHRD:

 Q = H(ilo) H(ilo+1) . . . H(ihi-1).

          N is INTEGER
          The order of the matrix Q. N >= 0.

          ILO is INTEGER

          IHI is INTEGER

          ILO and IHI must have the same values as in the previous call
          of SGEHRD. Q is equal to the unit matrix except in the
          submatrix Q(ilo+1:ihi,ilo+1:ihi).
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

          A is REAL array, dimension (LDA,N)
          On entry, the vectors which define the elementary reflectors,
          as returned by SGEHRD.
          On exit, the N-by-N orthogonal matrix Q.

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

          TAU is REAL array, dimension (N-1)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEHRD.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= IHI-ILO.
          For optimum performance LWORK >= (IHI-ILO)*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORGL2 generates an m by n real matrix Q with orthonormal rows,
 which is defined as the first m rows of a product of k elementary
 reflectors of order n

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

 as returned by SGELQF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. N >= M.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. M >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the i-th row must contain the vector which defines
          the elementary reflector H(i), for i = 1,2,...,k, as returned
          by SGELQF in the first k rows of its array argument A.
          On exit, the m-by-n matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGELQF.

          WORK is REAL array, dimension (M)

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument has an illegal value

 SORGLQ generates an M-by-N real matrix Q with orthonormal rows,
 which is defined as the first M rows of a product of K elementary
 reflectors of order N

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

 as returned by SGELQF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. N >= M.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. M >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the i-th row must contain the vector which defines
          the elementary reflector H(i), for i = 1,2,...,k, as returned
          by SGELQF in the first k rows of its array argument A.
          On exit, the M-by-N matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGELQF.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,M).
          For optimum performance LWORK >= M*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument has an illegal value

 SORGQL generates an M-by-N real matrix Q with orthonormal columns,
 which is defined as the last N columns of a product of K elementary
 reflectors of order M

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

 as returned by SGEQLF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. M >= N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. N >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the (n-k+i)-th column must contain the vector which
          defines the elementary reflector H(i), for i = 1,2,...,k, as
          returned by SGEQLF in the last k columns of its array
          argument A.
          On exit, the M-by-N matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQLF.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,N).
          For optimum performance LWORK >= N*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument has an illegal value

 SORGQR generates an M-by-N real matrix Q with orthonormal columns,
 which is defined as the first N columns of a product of K elementary
 reflectors of order M

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

 as returned by SGEQRF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. M >= N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. N >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the i-th column must contain the vector which
          defines the elementary reflector H(i), for i = 1,2,...,k, as
          returned by SGEQRF in the first k columns of its array
          argument A.
          On exit, the M-by-N matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQRF.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,N).
          For optimum performance LWORK >= N*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument has an illegal value

 SORGR2 generates an m by n real matrix Q with orthonormal rows,
 which is defined as the last m rows of a product of k elementary
 reflectors of order n

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

 as returned by SGERQF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. N >= M.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. M >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the (m-k+i)-th row must contain the vector which
          defines the elementary reflector H(i), for i = 1,2,...,k, as
          returned by SGERQF in the last k rows of its array argument
          A.
          On exit, the m by n matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGERQF.

          WORK is REAL array, dimension (M)

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument has an illegal value

 SORGRQ generates an M-by-N real matrix Q with orthonormal rows,
 which is defined as the last M rows of a product of K elementary
 reflectors of order N

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

 as returned by SGERQF.

          M is INTEGER
          The number of rows of the matrix Q. M >= 0.

          N is INTEGER
          The number of columns of the matrix Q. N >= M.

          K is INTEGER
          The number of elementary reflectors whose product defines the
          matrix Q. M >= K >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the (m-k+i)-th row must contain the vector which
          defines the elementary reflector H(i), for i = 1,2,...,k, as
          returned by SGERQF in the last k rows of its array argument
          A.
          On exit, the M-by-N matrix Q.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGERQF.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,M).
          For optimum performance LWORK >= M*NB, where NB is the
          optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument has an illegal value

 SORGTR generates a real orthogonal matrix Q which is defined as the
 product of n-1 elementary reflectors of order N, as returned by
 SSYTRD:

 if UPLO = 'U', Q = H(n-1) . . . H(2) H(1),

 if UPLO = 'L', Q = H(1) H(2) . . . H(n-1).

          UPLO is CHARACTER*1
          = 'U': Upper triangle of A contains elementary reflectors
                 from SSYTRD;
          = 'L': Lower triangle of A contains elementary reflectors
                 from SSYTRD.

          N is INTEGER
          The order of the matrix Q. N >= 0.

          A is REAL array, dimension (LDA,N)
          On entry, the vectors which define the elementary reflectors,
          as returned by SSYTRD.
          On exit, the N-by-N orthogonal matrix Q.

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

          TAU is REAL array, dimension (N-1)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SSYTRD.

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,N-1).
          For optimum performance LWORK >= (N-1)*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORM2L overwrites the general real m by n matrix C with

       Q * C  if SIDE = 'L' and TRANS = 'N', or

       Q**T * C  if SIDE = 'L' and TRANS = 'T', or

       C * Q  if SIDE = 'R' and TRANS = 'N', or

       C * Q**T if SIDE = 'R' and TRANS = 'T',

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGEQLF. Q is of order m if SIDE = 'L' and of order n
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left
          = 'R': apply Q or Q**T from the Right

          TRANS is CHARACTER*1
          = 'N': apply Q  (No transpose)
          = 'T': apply Q**T (Transpose)

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension (LDA,K)
          The i-th column must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGEQLF in the last k columns of its array argument A.
          A is modified by the routine but restored on exit.

          LDA is INTEGER
          The leading dimension of the array A.
          If SIDE = 'L', LDA >= max(1,M);
          if SIDE = 'R', LDA >= max(1,N).

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQLF.

          C is REAL array, dimension (LDC,N)
          On entry, the m by n matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                                   (N) if SIDE = 'L',
                                   (M) if SIDE = 'R'

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

 SORM2R overwrites the general real m by n matrix C with

       Q * C  if SIDE = 'L' and TRANS = 'N', or

       Q**T* C  if SIDE = 'L' and TRANS = 'T', or

       C * Q  if SIDE = 'R' and TRANS = 'N', or

       C * Q**T if SIDE = 'R' and TRANS = 'T',

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGEQRF. Q is of order m if SIDE = 'L' and of order n
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left
          = 'R': apply Q or Q**T from the Right

          TRANS is CHARACTER*1
          = 'N': apply Q  (No transpose)
          = 'T': apply Q**T (Transpose)

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension (LDA,K)
          The i-th column must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGEQRF in the first k columns of its array argument A.
          A is modified by the routine but restored on exit.

          LDA is INTEGER
          The leading dimension of the array A.
          If SIDE = 'L', LDA >= max(1,M);
          if SIDE = 'R', LDA >= max(1,N).

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQRF.

          C is REAL array, dimension (LDC,N)
          On entry, the m by n matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                                   (N) if SIDE = 'L',
                                   (M) if SIDE = 'R'

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

 If VECT = 'Q', SORMBR overwrites the general real M-by-N matrix C
 with
                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 If VECT = 'P', SORMBR overwrites the general real M-by-N matrix C
 with
                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      P * C          C * P
 TRANS = 'T':      P**T * C       C * P**T

 Here Q and P**T are the orthogonal matrices determined by SGEBRD when
 reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and
 P**T are defined as products of elementary reflectors H(i) and G(i)
 respectively.

 Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the
 order of the orthogonal matrix Q or P**T that is applied.

 If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:
 if nq >= k, Q = H(1) H(2) . . . H(k);
 if nq < k, Q = H(1) H(2) . . . H(nq-1).

 If VECT = 'P', A is assumed to have been a K-by-NQ matrix:
 if k < nq, P = G(1) G(2) . . . G(k);
 if k >= nq, P = G(1) G(2) . . . G(nq-1).

          VECT is CHARACTER*1
          = 'Q': apply Q or Q**T;
          = 'P': apply P or P**T.

          SIDE is CHARACTER*1
          = 'L': apply Q, Q**T, P or P**T from the Left;
          = 'R': apply Q, Q**T, P or P**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q  or P;
          = 'T':  Transpose, apply Q**T or P**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          If VECT = 'Q', the number of columns in the original
          matrix reduced by SGEBRD.
          If VECT = 'P', the number of rows in the original
          matrix reduced by SGEBRD.
          K >= 0.

          A is REAL array, dimension
                                (LDA,min(nq,K)) if VECT = 'Q'
                                (LDA,nq)        if VECT = 'P'
          The vectors which define the elementary reflectors H(i) and
          G(i), whose products determine the matrices Q and P, as
          returned by SGEBRD.

          LDA is INTEGER
          The leading dimension of the array A.
          If VECT = 'Q', LDA >= max(1,nq);
          if VECT = 'P', LDA >= max(1,min(nq,K)).

          TAU is REAL array, dimension (min(nq,K))
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i) or G(i) which determines Q or P, as returned
          by SGEBRD in the array argument TAUQ or TAUP.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q
          or P*C or P**T*C or C*P or C*P**T.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For optimum performance LWORK >= N*NB if SIDE = 'L', and
          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
          blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORMHR overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix of order nq, with nq = m if
 SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
 IHI-ILO elementary reflectors, as returned by SGEHRD:

 Q = H(ilo) H(ilo+1) . . . H(ihi-1).

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          ILO is INTEGER

          IHI is INTEGER

          ILO and IHI must have the same values as in the previous call
          of SGEHRD. Q is equal to the unit matrix except in the
          submatrix Q(ilo+1:ihi,ilo+1:ihi).
          If SIDE = 'L', then 1 <= ILO <= IHI <= M, if M > 0, and
          ILO = 1 and IHI = 0, if M = 0;
          if SIDE = 'R', then 1 <= ILO <= IHI <= N, if N > 0, and
          ILO = 1 and IHI = 0, if N = 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L'
                               (LDA,N) if SIDE = 'R'
          The vectors which define the elementary reflectors, as
          returned by SGEHRD.

          LDA is INTEGER
          The leading dimension of the array A.
          LDA >= max(1,M) if SIDE = 'L'; LDA >= max(1,N) if SIDE = 'R'.

          TAU is REAL array, dimension
                               (M-1) if SIDE = 'L'
                               (N-1) if SIDE = 'R'
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEHRD.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For optimum performance LWORK >= N*NB if SIDE = 'L', and
          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
          blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORML2 overwrites the general real m by n matrix C with

       Q * C  if SIDE = 'L' and TRANS = 'N', or

       Q**T* C  if SIDE = 'L' and TRANS = 'T', or

       C * Q  if SIDE = 'R' and TRANS = 'N', or

       C * Q**T if SIDE = 'R' and TRANS = 'T',

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGELQF. Q is of order m if SIDE = 'L' and of order n
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left
          = 'R': apply Q or Q**T from the Right

          TRANS is CHARACTER*1
          = 'N': apply Q  (No transpose)
          = 'T': apply Q**T (Transpose)

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L',
                               (LDA,N) if SIDE = 'R'
          The i-th row must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGELQF in the first k rows of its array argument A.
          A is modified by the routine but restored on exit.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGELQF.

          C is REAL array, dimension (LDC,N)
          On entry, the m by n matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                                   (N) if SIDE = 'L',
                                   (M) if SIDE = 'R'

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

 SORMLQ overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGELQF. Q is of order M if SIDE = 'L' and of order N
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L',
                               (LDA,N) if SIDE = 'R'
          The i-th row must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGELQF in the first k rows of its array argument A.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGELQF.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORMQL overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGEQLF. Q is of order M if SIDE = 'L' and of order N
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension (LDA,K)
          The i-th column must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGEQLF in the last k columns of its array argument A.

          LDA is INTEGER
          The leading dimension of the array A.
          If SIDE = 'L', LDA >= max(1,M);
          if SIDE = 'R', LDA >= max(1,N).

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQLF.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORMQR overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGEQRF. Q is of order M if SIDE = 'L' and of order N
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension (LDA,K)
          The i-th column must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGEQRF in the first k columns of its array argument A.

          LDA is INTEGER
          The leading dimension of the array A.
          If SIDE = 'L', LDA >= max(1,M);
          if SIDE = 'R', LDA >= max(1,N).

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGEQRF.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORMR2 overwrites the general real m by n matrix C with

       Q * C  if SIDE = 'L' and TRANS = 'N', or

       Q**T* C  if SIDE = 'L' and TRANS = 'T', or

       C * Q  if SIDE = 'R' and TRANS = 'N', or

       C * Q**T if SIDE = 'R' and TRANS = 'T',

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGERQF. Q is of order m if SIDE = 'L' and of order n
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left
          = 'R': apply Q or Q**T from the Right

          TRANS is CHARACTER*1
          = 'N': apply Q  (No transpose)
          = 'T': apply Q' (Transpose)

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L',
                               (LDA,N) if SIDE = 'R'
          The i-th row must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGERQF in the last k rows of its array argument A.
          A is modified by the routine but restored on exit.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGERQF.

          C is REAL array, dimension (LDC,N)
          On entry, the m by n matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                                   (N) if SIDE = 'L',
                                   (M) if SIDE = 'R'

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

 SORMR3 overwrites the general real m by n matrix C with

       Q * C  if SIDE = 'L' and TRANS = 'N', or

       Q**T* C  if SIDE = 'L' and TRANS = 'C', or

       C * Q  if SIDE = 'R' and TRANS = 'N', or

       C * Q**T if SIDE = 'R' and TRANS = 'C',

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by STZRZF. Q is of order m if SIDE = 'L' and of order n
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left
          = 'R': apply Q or Q**T from the Right

          TRANS is CHARACTER*1
          = 'N': apply Q  (No transpose)
          = 'T': apply Q**T (Transpose)

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          L is INTEGER
          The number of columns of the matrix A containing
          the meaningful part of the Householder reflectors.
          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L',
                               (LDA,N) if SIDE = 'R'
          The i-th row must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          STZRZF in the last k rows of its array argument A.
          A is modified by the routine but restored on exit.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by STZRZF.

          C is REAL array, dimension (LDC,N)
          On entry, the m-by-n matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension
                                   (N) if SIDE = 'L',
                                   (M) if SIDE = 'R'

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

 SORMRQ overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by SGERQF. Q is of order M if SIDE = 'L' and of order N
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L',
                               (LDA,N) if SIDE = 'R'
          The i-th row must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          SGERQF in the last k rows of its array argument A.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SGERQF.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORMRZ overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix defined as the product of k
 elementary reflectors

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

 as returned by STZRZF. Q is of order M if SIDE = 'L' and of order N
 if SIDE = 'R'.

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q.
          If SIDE = 'L', M >= K >= 0;
          if SIDE = 'R', N >= K >= 0.

          L is INTEGER
          The number of columns of the matrix A containing
          the meaningful part of the Householder reflectors.
          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L',
                               (LDA,N) if SIDE = 'R'
          The i-th row must contain the vector which defines the
          elementary reflector H(i), for i = 1,2,...,k, as returned by
          STZRZF in the last k rows of its array argument A.
          A is modified by the routine but restored on exit.

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

          TAU is REAL array, dimension (K)
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by STZRZF.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SORMTR overwrites the general real M-by-N matrix C with

                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T

 where Q is a real orthogonal matrix of order nq, with nq = m if
 SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
 nq-1 elementary reflectors, as returned by SSYTRD:

 if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);

 if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).

          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.

          UPLO is CHARACTER*1
          = 'U': Upper triangle of A contains elementary reflectors
                 from SSYTRD;
          = 'L': Lower triangle of A contains elementary reflectors
                 from SSYTRD.

          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.

          M is INTEGER
          The number of rows of the matrix C. M >= 0.

          N is INTEGER
          The number of columns of the matrix C. N >= 0.

          A is REAL array, dimension
                               (LDA,M) if SIDE = 'L'
                               (LDA,N) if SIDE = 'R'
          The vectors which define the elementary reflectors, as
          returned by SSYTRD.

          LDA is INTEGER
          The leading dimension of the array A.
          LDA >= max(1,M) if SIDE = 'L'; LDA >= max(1,N) if SIDE = 'R'.

          TAU is REAL array, dimension
                               (M-1) if SIDE = 'L'
                               (N-1) if SIDE = 'R'
          TAU(i) must contain the scalar factor of the elementary
          reflector H(i), as returned by SSYTRD.

          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).

          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.

          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For optimum performance LWORK >= N*NB if SIDE = 'L', and
          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
          blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

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

 SPBCON estimates the reciprocal of the condition number (in the
 1-norm) of a real symmetric positive definite band matrix using the
 Cholesky factorization A = U**T*U or A = L*L**T computed by SPBTRF.

 An estimate is obtained for norm(inv(A)), and the reciprocal of the
 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

          UPLO is CHARACTER*1
          = 'U':  Upper triangular factor stored in AB;
          = 'L':  Lower triangular factor stored in AB.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          AB is REAL array, dimension (LDAB,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T of the band matrix A, stored in the
          first KD+1 rows of the array.  The j-th column of U or L is
          stored in the j-th column of the array AB as follows:
          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).

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

          ANORM is REAL
          The 1-norm (or infinity-norm) of the symmetric band matrix A.

          RCOND is REAL
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
          estimate of the 1-norm of inv(A) computed in this routine.

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

 SPBEQU computes row and column scalings intended to equilibrate a
 symmetric positive definite band matrix A and reduce its condition
 number (with respect to the two-norm).  S contains the scale factors,
 S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
 elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal.  This
 choice of S puts the condition number of B within a factor N of the
 smallest possible condition number over all possible diagonal
 scalings.

          UPLO is CHARACTER*1
          = 'U':  Upper triangular of A is stored;
          = 'L':  Lower triangular of A is stored.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          AB is REAL array, dimension (LDAB,N)
          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).

          LDAB is INTEGER
          The leading dimension of the array A.  LDAB >= KD+1.

          S is REAL array, dimension (N)
          If INFO = 0, S contains the scale factors for A.

          SCOND is REAL
          If INFO = 0, S contains the ratio of the smallest S(i) to
          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
          large nor too small, it is not worth scaling by S.

          AMAX is REAL
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if INFO = i, the i-th diagonal element is nonpositive.

 SPBRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite
 and banded, and provides error bounds and backward error estimates
 for the solution.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.

          AB is REAL array, dimension (LDAB,N)
          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).

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

          AFB is REAL array, dimension (LDAFB,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T of the band matrix A as computed by
          SPBTRF, in the same storage format as A (see AB).

          LDAFB is INTEGER
          The leading dimension of the array AFB.  LDAFB >= KD+1.

          B is REAL array, dimension (LDB,NRHS)
          The right hand side matrix B.

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

          X is REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SPBTRS.
          On exit, the improved solution matrix X.

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

          FERR is REAL array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.

          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

  ITMAX is the maximum number of steps of iterative refinement.

 SPBSTF computes a split Cholesky factorization of a real
 symmetric positive definite band matrix A.

 This routine is designed to be used in conjunction with SSBGST.

 The factorization has the form  A = S**T*S  where S is a band matrix
 of the same bandwidth as A and the following structure:

   S = ( U    )
       ( M  L )

 where U is upper triangular of order m = (n+kd)/2, and L is lower
 triangular of order n-m.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          AB is REAL 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 factor S from the split Cholesky
          factorization A = S**T*S. See Further Details.

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

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, the factorization could not be completed,
               because the updated element a(i,i) was negative; the
               matrix A is not positive definite.

  The band storage scheme is illustrated by the following example, when
  N = 7, KD = 2:

  S = ( s11  s12  s13                     )
      (      s22  s23  s24                )
      (           s33  s34                )
      (                s44                )
      (           s53  s54  s55           )
      (                s64  s65  s66      )
      (                     s75  s76  s77 )

  If UPLO = 'U', the array AB holds:

  on entry:                          on exit:

   *    *   a13  a24  a35  a46  a57   *    *   s13  s24  s53  s64  s75
   *   a12  a23  a34  a45  a56  a67   *   s12  s23  s34  s54  s65  s76
  a11  a22  a33  a44  a55  a66  a77  s11  s22  s33  s44  s55  s66  s77

  If UPLO = 'L', the array AB holds:

  on entry:                          on exit:

  a11  a22  a33  a44  a55  a66  a77  s11  s22  s33  s44  s55  s66  s77
  a21  a32  a43  a54  a65  a76   *   s12  s23  s34  s54  s65  s76   *
  a31  a42  a53  a64  a64   *    *   s13  s24  s53  s64  s75   *    *

  Array elements marked * are not used by the routine.

 SPBTF2 computes the Cholesky factorization of a real symmetric
 positive definite band matrix A.

 The factorization has the form
    A = U**T * U ,  if UPLO = 'U', or
    A = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix, U**T is the transpose of U, and
 L is lower triangular.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.

          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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 is REAL 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 is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -k, the k-th argument had an illegal value
          > 0: if INFO = k, the leading minor of order k is not
               positive definite, and the factorization could not be
               completed.

  The band storage scheme is illustrated by the following example, when
  N = 6, KD = 2, and UPLO = 'U':

  On entry:                       On exit:

      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66

  Similarly, if UPLO = 'L' the format of A is as follows:

  On entry:                       On exit:

     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *

  Array elements marked * are not used by the routine.

 SPBTRF computes the Cholesky factorization of a real symmetric
 positive definite band matrix A.

 The factorization has the form
    A = U**T * U,  if UPLO = 'U', or
    A = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          AB is REAL 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 is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the leading minor of order i is not
                positive definite, and the factorization could not be
                completed.

  The band storage scheme is illustrated by the following example, when
  N = 6, KD = 2, and UPLO = 'U':

  On entry:                       On exit:

      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66

  Similarly, if UPLO = 'L' the format of A is as follows:

  On entry:                       On exit:

     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *

  Array elements marked * are not used by the routine.

 SPBTRS solves a system of linear equations A*X = B with a symmetric
 positive definite band matrix A using the Cholesky factorization
 A = U**T*U or A = L*L**T computed by SPBTRF.

          UPLO is CHARACTER*1
          = 'U':  Upper triangular factor stored in AB;
          = 'L':  Lower triangular factor stored in AB.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.

          AB is REAL array, dimension (LDAB,N)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T of the band matrix A, stored in the
          first KD+1 rows of the array.  The j-th column of U or L is
          stored in the j-th column of the array AB as follows:
          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).

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

          B is REAL array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.

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

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

 SPFTRF computes the Cholesky factorization of a real symmetric
 positive definite matrix A.

 The factorization has the form
    A = U**T * U,  if UPLO = 'U', or
    A = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular.

 This is the block version of the algorithm, calling Level 3 BLAS.

          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of RFP A is stored;
          = 'L':  Lower triangle of RFP A is stored.

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

          A is REAL array, dimension ( N*(N+1)/2 );
          On entry, the symmetric matrix A in RFP format. RFP format is
          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
          the transpose of RFP A as defined when
          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
          follows: If UPLO = 'U' the RFP A contains the NT elements of
          upper packed A. If UPLO = 'L' the RFP A contains the elements
          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
          'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
          is odd. See the Note below for more details.

          On exit, if INFO = 0, the factor U or L from the Cholesky
          factorization RFP A = U**T*U or RFP A = L*L**T.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the leading minor of order i is not
                positive definite, and the factorization could not be
                completed.

  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

 SPFTRI computes the inverse of a real (symmetric) positive definite
 matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
 computed by SPFTRF.

          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          A is REAL array, dimension ( N*(N+1)/2 )
          On entry, the symmetric matrix A in RFP format. RFP format is
          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
          the transpose of RFP A as defined when
          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
          follows: If UPLO = 'U' the RFP A contains the nt elements of
          upper packed A. If UPLO = 'L' the RFP A contains the elements
          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
          'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
          is odd. See the Note below for more details.

          On exit, the symmetric inverse of the original matrix, in the
          same storage format.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the (i,i) element of the factor U or L is
                zero, and the inverse could not be computed.

  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

 SPFTRS solves a system of linear equations A*X = B with a symmetric
 positive definite matrix A using the Cholesky factorization
 A = U**T*U or A = L*L**T computed by SPFTRF.

          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of RFP A is stored;
          = 'L':  Lower triangle of RFP A is stored.

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

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.

          A is REAL array, dimension ( N*(N+1)/2 )
          The triangular factor U or L from the Cholesky factorization
          of RFP A = U**H*U or RFP A = L*L**T, as computed by SPFTRF.
          See note below for more details about RFP A.

          B is REAL array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.

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

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

  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

 SPPCON estimates the reciprocal of the condition number (in the
 1-norm) of a real symmetric positive definite packed matrix using
 the Cholesky factorization A = U**T*U or A = L*L**T computed by
 SPPTRF.

 An estimate is obtained for norm(inv(A)), and the reciprocal of the
 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          AP is REAL array, dimension (N*(N+1)/2)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T, packed columnwise in a linear
          array.  The j-th column of U or L is stored in the array AP
          as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.

          ANORM is REAL
          The 1-norm (or infinity-norm) of the symmetric matrix A.

          RCOND is REAL
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
          estimate of the 1-norm of inv(A) computed in this routine.

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

 SPPEQU computes row and column scalings intended to equilibrate a
 symmetric positive definite matrix A in packed storage and reduce
 its condition number (with respect to the two-norm).  S contains the
 scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix
 B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal.
 This choice of S puts the condition number of B within a factor N of
 the smallest possible condition number over all possible diagonal
 scalings.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          AP is REAL 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.

          S is REAL array, dimension (N)
          If INFO = 0, S contains the scale factors for A.

          SCOND is REAL
          If INFO = 0, S contains the ratio of the smallest S(i) to
          the largest S(i).  If SCOND >= 0.1 and AMAX is neither too
          large nor too small, it is not worth scaling by S.

          AMAX is REAL
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the i-th diagonal element is nonpositive.

 SPPRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric positive definite
 and packed, and provides error bounds and backward error estimates
 for the solution.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.

          AP is REAL 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.

          AFP is REAL array, dimension (N*(N+1)/2)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T, as computed by SPPTRF/CPPTRF,
          packed columnwise in a linear array in the same format as A
          (see AP).

          B is REAL array, dimension (LDB,NRHS)
          The right hand side matrix B.

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

          X is REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SPPTRS.
          On exit, the improved solution matrix X.

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

          FERR is REAL array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.

          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

  ITMAX is the maximum number of steps of iterative refinement.

 SPPTRF computes the Cholesky factorization of a real symmetric
 positive definite matrix A stored in packed format.

 The factorization has the form
    A = U**T * U,  if UPLO = 'U', or
    A = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          AP is REAL 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.
          See below for further details.

          On exit, if INFO = 0, the triangular factor U or L from the
          Cholesky factorization A = U**T*U or A = L*L**T, in the same
          storage format as A.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the leading minor of order i is not
                positive definite, and the factorization could not be
                completed.

  The packed storage scheme is illustrated by the following example
  when N = 4, UPLO = 'U':

  Two-dimensional storage of the symmetric matrix A:

     a11 a12 a13 a14
         a22 a23 a24
             a33 a34     (aij = aji)
                 a44

  Packed storage of the upper triangle of A:

  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]

 SPPTRI computes the inverse of a real symmetric positive definite
 matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
 computed by SPPTRF.

          UPLO is CHARACTER*1
          = 'U':  Upper triangular factor is stored in AP;
          = 'L':  Lower triangular factor is stored in AP.

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

          AP is REAL array, dimension (N*(N+1)/2)
          On entry, the triangular factor U or L from the Cholesky
          factorization A = U**T*U or A = L*L**T, packed columnwise as
          a linear array.  The j-th column of U or L is stored in the
          array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.

          On exit, the upper or lower triangle of the (symmetric)
          inverse of A, overwriting the input factor U or L.

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the (i,i) element of the factor U or L is
                zero, and the inverse could not be computed.

 SPPTRS solves a system of linear equations A*X = B with a symmetric
 positive definite matrix A in packed storage using the Cholesky
 factorization A = U**T*U or A = L*L**T computed by SPPTRF.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.

          AP is REAL array, dimension (N*(N+1)/2)
          The triangular factor U or L from the Cholesky factorization
          A = U**T*U or A = L*L**T, packed columnwise in a linear
          array.  The j-th column of U or L is stored in the array AP
          as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.

          B is REAL array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.

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

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

 SPSTF2 computes the Cholesky factorization with complete
 pivoting of a real symmetric positive semidefinite matrix A.

 The factorization has the form
    P**T * A * P = U**T * U ,  if UPLO = 'U',
    P**T * A * P = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular, and
 P is stored as vector PIV.

 This algorithm does not attempt to check that A is positive
 semidefinite. This version of the algorithm calls level 2 BLAS.

          UPLO 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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 INFO = 0, the factor U or L from the Cholesky
          factorization as above.

          PIV is INTEGER array, dimension (N)
          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.

          RANK is INTEGER
          The rank of A given by the number of steps the algorithm
          completed.

          TOL is REAL
          User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
          will be used. The algorithm terminates at the (K-1)st step
          if the pivot <= TOL.

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

          WORK is REAL array, dimension (2*N)
          Work space.

          INFO is INTEGER
          < 0: If INFO = -K, the K-th argument had an illegal value,
          = 0: algorithm completed successfully, and
          > 0: the matrix A is either rank deficient with computed rank
               as returned in RANK, or is not positive semidefinite. See
               Section 7 of LAPACK Working Note #161 for further
               information.

 SPSTRF computes the Cholesky factorization with complete
 pivoting of a real symmetric positive semidefinite matrix A.

 The factorization has the form
    P**T * A * P = U**T * U ,  if UPLO = 'U',
    P**T * A * P = L  * L**T,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular, and
 P is stored as vector PIV.

 This algorithm does not attempt to check that A is positive
 semidefinite. This version of the algorithm calls level 3 BLAS.

          UPLO 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 is INTEGER
          The order of the matrix A.  N >= 0.

          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 INFO = 0, the factor U or L from the Cholesky
          factorization as above.

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

          PIV is INTEGER array, dimension (N)
          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.

          RANK is INTEGER
          The rank of A given by the number of steps the algorithm
          completed.

          TOL is REAL
          User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
          will be used. The algorithm terminates at the (K-1)st step
          if the pivot <= TOL.

          WORK is REAL array, dimension (2*N)
          Work space.

          INFO is INTEGER
          < 0: If INFO = -K, the K-th argument had an illegal value,
          = 0: algorithm completed successfully, and
          > 0: the matrix A is either rank deficient with computed rank
               as returned in RANK, or is not positive semidefinite. See
               Section 7 of LAPACK Working Note #161 for further
               information.

 SSBGST reduces a real symmetric-definite banded generalized
 eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
 such that C has the same bandwidth as A.

 B must have been previously factorized as S**T*S by SPBSTF, using a
 split Cholesky factorization. A is overwritten by C = X**T*A*X, where
 X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
 bandwidth of A.

          VECT is CHARACTER*1
          = 'N':  do not form the transformation matrix X;
          = 'V':  form X.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          KA is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.

          KB is INTEGER
          The number of superdiagonals of the matrix B if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.

          AB is REAL array, dimension (LDAB,N)
          On entry, the upper or lower triangle of the symmetric band
          matrix A, stored in the first ka+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).

          On exit, the transformed matrix X**T*A*X, stored in the same
          format as A.

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

          BB is REAL array, dimension (LDBB,N)
          The banded factor S from the split Cholesky factorization of
          B, as returned by SPBSTF, stored in the first KB+1 rows of
          the array.

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

          X is REAL array, dimension (LDX,N)
          If VECT = 'V', the n-by-n matrix X.
          If VECT = 'N', the array X is not referenced.

          LDX is INTEGER
          The leading dimension of the array X.
          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.

          WORK is REAL array, dimension (2*N)

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

 SSBTRD reduces a real symmetric band matrix A to symmetric
 tridiagonal form T by an orthogonal similarity transformation:
 Q**T * A * Q = T.

          VECT is CHARACTER*1
          = 'N':  do not form Q;
          = 'V':  form Q;
          = 'U':  update a matrix X, by forming X*Q.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.

          AB is REAL 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, the diagonal elements of AB are overwritten by the
          diagonal elements of the tridiagonal matrix T; if KD > 0, the
          elements on the first superdiagonal (if UPLO = 'U') or the
          first subdiagonal (if UPLO = 'L') are overwritten by the
          off-diagonal elements of T; the rest of AB is overwritten by
          values generated during the reduction.

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

          D is REAL array, dimension (N)
          The diagonal elements of the tridiagonal matrix T.

          E is REAL array, dimension (N-1)
          The off-diagonal elements of the tridiagonal matrix T:
          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.

          Q is REAL array, dimension (LDQ,N)
          On entry, if VECT = 'U', then Q must contain an N-by-N
          matrix X; if VECT = 'N' or 'V', then Q need not be set.

          On exit:
          if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
          if VECT = 'U', Q contains the product X*Q;
          if VECT = 'N', the array Q is not referenced.

          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.

          WORK is REAL array, dimension (N)

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

  Modified by Linda Kaufman, Bell Labs.

 Level 3 BLAS like routine for C in RFP Format.

 SSFRK performs one of the symmetric rank--k operations

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

 or

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

 where alpha and beta are real scalars, C is an n--by--n symmetric
 matrix and A is an n--by--k matrix in the first case and a k--by--n
 matrix in the second case.

          TRANSR is CHARACTER*1
          = 'N':  The Normal Form of RFP A is stored;
          = 'T':  The Transpose Form of RFP A is stored.

          UPLO is CHARACTER*1
           On  entry, UPLO specifies whether the upper or lower
           triangular part of the array C is to be referenced as
           follows:

              UPLO = 'U' or 'u'   Only the upper triangular part of C
                                  is to be referenced.

              UPLO = 'L' or 'l'   Only the lower triangular part of C
                                  is to be referenced.

           Unchanged on exit.

          TRANS is CHARACTER*1
           On entry, TRANS specifies the operation to be performed as
           follows:

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

              TRANS = 'T' or 't'   C := alpha*A**T*A + beta*C.

           Unchanged on exit.

          N is INTEGER
           On entry, N specifies the order of the matrix C. N must be
           at least zero.
           Unchanged on exit.

          K is INTEGER
           On entry with TRANS = 'N' or 'n', K specifies the number
           of  columns of the matrix A, and on entry with TRANS = 'T'
           or 't', K specifies the number of rows of the matrix A. K
           must be at least zero.
           Unchanged on exit.

          ALPHA is REAL
           On entry, ALPHA specifies the scalar alpha.
           Unchanged on exit.

          A is REAL array of DIMENSION (LDA,ka)
           where KA
           is K  when TRANS = 'N' or 'n', and is N otherwise. Before
           entry with TRANS = 'N' or 'n', the leading N--by--K part of
           the array A must contain the matrix A, otherwise the leading
           K--by--N part of the array A must contain the matrix A.
           Unchanged on exit.

          LDA is INTEGER
           On entry, LDA specifies the first dimension of A as declared
           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
           then  LDA must be at least  max( 1, n ), otherwise  LDA must
           be at least  max( 1, k ).
           Unchanged on exit.

          BETA is REAL
           On entry, BETA specifies the scalar beta.
           Unchanged on exit.

          C is REAL array, dimension (NT)
           NT = N*(N+1)/2. On entry, the symmetric matrix C in RFP
           Format. RFP Format is described by TRANSR, UPLO and N.

 SSPCON estimates the reciprocal of the condition number (in the
 1-norm) of a real symmetric packed matrix A using the factorization
 A = U*D*U**T or A = L*D*L**T computed by SSPTRF.

 An estimate is obtained for norm(inv(A)), and the reciprocal of the
 condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).

          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.

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

          AP is REAL array, dimension (N*(N+1)/2)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by SSPTRF, stored as a
          packed triangular matrix.

          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by SSPTRF.

          ANORM is REAL
          The 1-norm of the original matrix A.

          RCOND is REAL
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
          estimate of the 1-norm of inv(A) computed in this routine.

          WORK is REAL array, dimension (2*N)

          IWORK is INTEGER array, dimension (N)

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

 SSPGST reduces a real symmetric-definite generalized eigenproblem
 to standard form, using packed storage.

 If ITYPE = 1, the problem is A*x = lambda*B*x,
 and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)

 If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
 B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.

 B must have been previously factorized as U**T*U or L*L**T by SPPTRF.

          ITYPE is INTEGER
          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
          = 2 or 3: compute U*A*U**T or L**T*A*L.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored and B is factored as
                  U**T*U;
          = 'L':  Lower triangle of A is stored and B is factored as
                  L*L**T.

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

          AP is REAL 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, if INFO = 0, the transformed matrix, stored in the
          same format as A.

          BP is REAL array, dimension (N*(N+1)/2)
          The triangular factor from the Cholesky factorization of B,
          stored in the same format as A, as returned by SPPTRF.

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

 SSPRFS improves the computed solution to a system of linear
 equations when the coefficient matrix is symmetric indefinite
 and packed, and provides error bounds and backward error estimates
 for the solution.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.

          AP is REAL array, dimension (N*(N+1)/2)
          The upper or lower triangle of the symmetric matrix A, packed
          columnwise in a linear array.  The j-th column of A is stored
          in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

          AFP is REAL array, dimension (N*(N+1)/2)
          The factored form of the matrix A.  AFP contains the block
          diagonal matrix D and the multipliers used to obtain the
          factor U or L from the factorization A = U*D*U**T or
          A = L*D*L**T as computed by SSPTRF, stored as a packed
          triangular matrix.

          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by SSPTRF.

          B is REAL array, dimension (LDB,NRHS)
          The right hand side matrix B.

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

          X is REAL array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by SSPTRS.
          On exit, the improved solution matrix X.

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

          FERR is REAL array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.

          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

  ITMAX is the maximum number of steps of iterative refinement.

 SSPTRD reduces a real symmetric matrix A stored in packed form to
 symmetric tridiagonal form T by an orthogonal similarity
 transformation: Q**T * A * Q = T.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          AP is REAL 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)*(2*n-j)/2) = A(i,j) for j<=i<=n.
          On exit, if UPLO = 'U', the diagonal and first superdiagonal
          of A are overwritten by the corresponding elements of the
          tridiagonal matrix T, and the elements above the first
          superdiagonal, with the array TAU, represent the orthogonal
          matrix Q as a product of elementary reflectors; if UPLO
          = 'L', the diagonal and first subdiagonal of A are over-
          written by the corresponding elements of the tridiagonal
          matrix T, and the elements below the first subdiagonal, with
          the array TAU, represent the orthogonal matrix Q as a product
          of elementary reflectors. See Further Details.

          D is REAL array, dimension (N)
          The diagonal elements of the tridiagonal matrix T:
          D(i) = A(i,i).

          E is REAL array, dimension (N-1)
          The off-diagonal elements of the tridiagonal matrix T:
          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.

          TAU is REAL array, dimension (N-1)
          The scalar factors of the elementary reflectors (see Further
          Details).

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

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

     Q = H(n-1) . . . H(2) H(1).

  Each H(i) has the form

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

  where tau is a real scalar, and v is a real vector with
  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).

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

     Q = H(1) H(2) . . . H(n-1).

  Each H(i) has the form

     H(i) = I - tau * v * v**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+2:n) is stored on exit in AP,
  overwriting A(i+2:n,i), and tau is stored in TAU(i).

 SSPTRF computes the factorization of a real symmetric matrix A stored
 in packed format using the Bunch-Kaufman diagonal pivoting method:

    A = U*D*U**T  or  A = L*D*L**T

 where U (or L) is a product of permutation and unit upper (lower)
 triangular matrices, and D is symmetric and block diagonal with
 1-by-1 and 2-by-2 diagonal blocks.

          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.

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

          AP is REAL array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the symmetric matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.

          On exit, the block diagonal matrix D and the multipliers used
          to obtain the factor U or L, stored as a packed triangular
          matrix overwriting A (see below for further details).

          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D.
          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
          interchanged and D(k,k) is a 1-by-1 diagonal block.
          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
               has been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if it
               is used to solve a system of equations.

  5-96 - Based on modifications by J. Lewis, Boeing Computer Services
         Company

  If UPLO = 'U', then A = U*D*U**T, where
     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

             (   I    v    0   )   k-s
     U(k) =  (   0    I    0   )   s
             (   0    0    I   )   n-k
                k-s   s   n-k

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

  If UPLO = 'L', then A = L*D*L**T, where
     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  that if the diagonal block D(k) is of order s (s = 1 or 2), then

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

  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).

 SSPTRI computes the inverse of a real symmetric indefinite matrix
 A in packed storage using the factorization A = U*D*U**T or
 A = L*D*L**T computed by SSPTRF.

          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.

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

          AP is REAL array, dimension (N*(N+1)/2)
          On entry, the block diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by SSPTRF,
          stored as a packed triangular matrix.

          On exit, if INFO = 0, the (symmetric) inverse of the original
          matrix, stored as a packed triangular matrix. The j-th column
          of inv(A) is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
          if UPLO = 'L',
             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.

          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by SSPTRF.

          WORK is REAL array, dimension (N)

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
               inverse could not be computed.

 SSPTRS solves a system of linear equations A*X = B with a real
 symmetric matrix A stored in packed format using the factorization
 A = U*D*U**T or A = L*D*L**T computed by SSPTRF.

          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.

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

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.

          AP is REAL array, dimension (N*(N+1)/2)
          The block diagonal matrix D and the multipliers used to
          obtain the factor U or L as computed by SSPTRF, stored as a
          packed triangular matrix.

          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the block structure of D
          as determined by SSPTRF.

          B is REAL array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, the solution matrix X.

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

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

 SSTEGR computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
 a well defined set of pairwise different real eigenvalues, the corresponding
 real eigenvectors are pairwise orthogonal.

 The spectrum may be computed either completely or partially by specifying
 either an interval (VL,VU] or a range of indices IL:IU for the desired
 eigenvalues.

 SSTEGR is a compatibility wrapper around the improved SSTEMR routine.
 See SSTEMR for further details.

 One important change is that the ABSTOL parameter no longer provides any
 benefit and hence is no longer used.

 Note : SSTEGR and SSTEMR work only on machines which follow
 IEEE-754 floating-point standard in their handling of infinities and
 NaNs.  Normal execution may create these exceptiona values and hence
 may abort due to a floating point exception in environments which
 do not conform to the IEEE-754 standard.

          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.

          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found.
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th eigenvalues will be found.

          N is INTEGER
          The order of the matrix.  N >= 0.

          D is REAL array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal matrix
          T. On exit, D is overwritten.

          E is REAL array, dimension (N)
          On entry, the (N-1) subdiagonal elements of the tridiagonal
          matrix T in elements 1 to N-1 of E. E(N) need not be set on
          input, but is used internally as workspace.
          On exit, E is overwritten.

          VL is REAL

          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.

          VU is REAL

          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.

          IL is INTEGER

          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0.
          Not referenced if RANGE = 'A' or 'V'.

          IU is INTEGER

          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0.
          Not referenced if RANGE = 'A' or 'V'.

          ABSTOL is REAL
          Unused.  Was the absolute error tolerance for the
          eigenvalues/eigenvectors in previous versions.

          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

          W is REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.

          Z is REAL array, dimension (LDZ, max(1,M) )
          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
          contain the orthonormal eigenvectors of the matrix T
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          If JOBZ = 'N', then Z is not referenced.
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and an upper bound must be used.
          Supplying N columns is always safe.

          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', then LDZ >= max(1,N).

          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 computed eigenvector
          is nonzero only in elements ISUPPZ( 2*i-1 ) through
          ISUPPZ( 2*i ). This is relevant in the case when the matrix
          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.

          WORK is REAL array, dimension (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal
          (and minimal) LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,18*N)
          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          IWORK is INTEGER array, dimension (LIWORK)
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

          LIWORK is INTEGER
          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
          if the eigenvectors are desired, and LIWORK >= max(1,8*N)
          if only the eigenvalues are to be computed.
          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal size of the IWORK array,
          returns this value as the first entry of the IWORK array, and
          no error message related to LIWORK is issued by XERBLA.

          INFO is INTEGER
          On exit, INFO
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = 1X, internal error in SLARRE,
                if INFO = 2X, internal error in SLARRV.
                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
                the nonzero error code returned by SLARRE or
                SLARRV, respectively.

 SSTEIN computes the eigenvectors of a real symmetric tridiagonal
 matrix T corresponding to specified eigenvalues, using inverse
 iteration.

 The maximum number of iterations allowed for each eigenvector is
 specified by an internal parameter MAXITS (currently set to 5).

          N is INTEGER
          The order of the matrix.  N >= 0.

          D is REAL array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.

          E is REAL array, dimension (N-1)
          The (n-1) subdiagonal elements of the tridiagonal matrix
          T, in elements 1 to N-1.

          M is INTEGER
          The number of eigenvectors to be found.  0 <= M <= N.

          W is REAL array, dimension (N)
          The first M elements of W contain the 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 SSTEBZ with ORDER = 'B' is expected here. )

          IBLOCK is INTEGER array, dimension (N)
          The submatrix indices associated with the corresponding
          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
          the first submatrix from the top, =2 if W(i) belongs to
          the second submatrix, etc.  ( The output array IBLOCK
          from SSTEBZ is expected here. )

          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into submatrices.
          The first submatrix consists of rows/columns 1 to
          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
          through ISPLIT( 2 ), etc.
          ( The output array ISPLIT from SSTEBZ is expected here. )

          Z is REAL array, dimension (LDZ, M)
          The computed eigenvectors.  The eigenvector associated
          with the eigenvalue W(i) is stored in the i-th column of
          Z.  Any vector which fails to converge is set to its current
          iterate after MAXITS iterations.

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

          WORK is REAL array, dimension (5*N)

          IWORK is INTEGER array, dimension (N)

          IFAIL is INTEGER array, dimension (M)
          On normal exit, all elements of IFAIL are zero.
          If one or more eigenvectors fail to converge after
          MAXITS iterations, then their indices are stored in
          array IFAIL.

          INFO is INTEGER
          = 0: successful exit.
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, then i eigenvectors failed to converge
               in MAXITS iterations.  Their indices are stored in
               array IFAIL.

  MAXITS  INTEGER, default = 5
          The maximum number of iterations performed.

  EXTRA   INTEGER, default = 2
          The number of iterations performed after norm growth
          criterion is satisfied, should be at least 1.

 SSTEMR computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
 a well defined set of pairwise different real eigenvalues, the corresponding
 real eigenvectors are pairwise orthogonal.

 The spectrum may be computed either completely or partially by specifying
 either an interval (VL,VU] or a range of indices IL:IU for the desired
 eigenvalues.

 Depending on the number of desired eigenvalues, these are computed either
 by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
 computed by the use of various suitable L D L^T factorizations near clusters
 of close eigenvalues (referred to as RRRs, Relatively Robust
 Representations). An informal sketch of the algorithm follows.

 For each unreduced block (submatrix) of T,
    (a) Compute T - sigma I  = L D L^T, so that L and D
        define all the wanted eigenvalues to high relative accuracy.
        This means that small relative changes in the entries of D and L
        cause only small relative changes in the eigenvalues and
        eigenvectors. The standard (unfactored) representation of the
        tridiagonal matrix T does not have this property in general.
    (b) Compute the eigenvalues to suitable accuracy.
        If the eigenvectors are desired, the algorithm attains full
        accuracy of the computed eigenvalues only right before
        the corresponding vectors have to be computed, see steps c) and d).
    (c) For each cluster of close eigenvalues, select a new
        shift close to the cluster, find a new factorization, and refine
        the shifted eigenvalues to suitable accuracy.
    (d) For each eigenvalue with a large enough relative separation compute
        the corresponding eigenvector by forming a rank revealing twisted
        factorization. Go back to (c) for any clusters that remain.

 For more details, see:
 - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
   to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
 - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
   Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
   2004.  Also LAPACK Working Note 154.
 - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
   tridiagonal eigenvalue/eigenvector problem",
   Computer Science Division Technical Report No. UCB/CSD-97-971,
   UC Berkeley, May 1997.

 Further Details
 1.SSTEMR works only on machines which follow IEEE-754
 floating-point standard in their handling of infinities and NaNs.
 This permits the use of efficient inner loops avoiding a check for
 zero divisors.

          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.

          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found.
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th eigenvalues will be found.

          N is INTEGER
          The order of the matrix.  N >= 0.

          D is REAL array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal matrix
          T. On exit, D is overwritten.

          E is REAL array, dimension (N)
          On entry, the (N-1) subdiagonal elements of the tridiagonal
          matrix T in elements 1 to N-1 of E. E(N) need not be set on
          input, but is used internally as workspace.
          On exit, E is overwritten.

          VL is REAL

          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.

          VU is REAL

          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.

          IL is INTEGER

          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0.
          Not referenced if RANGE = 'A' or 'V'.

          IU is INTEGER

          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0.
          Not referenced if RANGE = 'A' or 'V'.

          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.

          W is REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.

          Z is REAL array, dimension (LDZ, max(1,M) )
          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
          contain the orthonormal eigenvectors of the matrix T
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          If JOBZ = 'N', then Z is not referenced.
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and can be computed with a workspace
          query by setting NZC = -1, see below.

          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', then LDZ >= max(1,N).

          NZC is INTEGER
          The number of eigenvectors to be held in the array Z.
          If RANGE = 'A', then NZC >= max(1,N).
          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
          If RANGE = 'I', then NZC >= IU-IL+1.
          If NZC = -1, then a workspace query is assumed; the
          routine calculates the number of columns of the array Z that
          are needed to hold the eigenvectors.
          This value is returned as the first entry of the Z array, and
          no error message related to NZC is issued by XERBLA.

          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 computed eigenvector
          is nonzero only in elements ISUPPZ( 2*i-1 ) through
          ISUPPZ( 2*i ). This is relevant in the case when the matrix
          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.

          TRYRAC is LOGICAL
          If TRYRAC.EQ..TRUE., indicates that the code should check whether
          the tridiagonal matrix defines its eigenvalues to high relative
          accuracy.  If so, the code uses relative-accuracy preserving
          algorithms that might be (a bit) slower depending on the matrix.
          If the matrix does not define its eigenvalues to high relative
          accuracy, the code can uses possibly faster algorithms.
          If TRYRAC.EQ..FALSE., the code is not required to guarantee
          relatively accurate eigenvalues and can use the fastest possible
          techniques.
          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
          does not define its eigenvalues to high relative accuracy.

          WORK is REAL array, dimension (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal
          (and minimal) LWORK.

          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,18*N)
          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.

          IWORK is INTEGER array, dimension (LIWORK)
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.

          LIWORK is INTEGER
          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
          if the eigenvectors are desired, and LIWORK >= max(1,8*N)
          if only the eigenvalues are to be computed.
          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal size of the IWORK array,
          returns this value as the first entry of the IWORK array, and
          no error message related to LIWORK is issued by XERBLA.

          INFO is INTEGER
          On exit, INFO
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = 1X, internal error in SLARRE,
                if INFO = 2X, internal error in SLARRV.
                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
                the nonzero error code returned by SLARRE or
                SLARRV, respectively.

 STBCON estimates the reciprocal of the condition number of a
 triangular band matrix A, in either the 1-norm or the infinity-norm.

 The norm of A is computed and an estimate is obtained for
 norm(inv(A)), then the reciprocal of the condition number is
 computed as
    RCOND = 1 / ( norm(A) * norm(inv(A)) ).

          NORM is CHARACTER*1
          Specifies whether the 1-norm condition number or the
          infinity-norm condition number is required:
          = '1' or 'O':  1-norm;
          = 'I':         Infinity-norm.

          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.

          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.

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

          KD is INTEGER
          The number of superdiagonals or subdiagonals of the
          triangular band matrix A.  KD >= 0.

          AB is REAL 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).
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.

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

          RCOND is REAL
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(norm(A) * norm(inv(A))).

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

 STBRFS provides error bounds and backward error estimates for the
 solution to a system of linear equations with a triangular band
 coefficient matrix.

 The solution matrix X must be computed by STBTRS or some other
 means before entering this routine.  STBRFS does not do iterative
 refinement because doing so cannot improve the backward error.

          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.

          TRANS is CHARACTER*1
          Specifies the form of the system of equations:
          = 'N':  A * X = B  (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)

          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.

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

          KD is INTEGER
          The number of superdiagonals or subdiagonals of the
          triangular band matrix A.  KD >= 0.

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X.  NRHS >= 0.

          AB is REAL 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).
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.

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

          B is REAL array, dimension (LDB,NRHS)
          The right hand side matrix B.

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

          X is REAL array, dimension (LDX,NRHS)
          The solution matrix X.

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

          FERR is REAL array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.

          BERR is REAL array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).

          WORK is REAL array, dimension (3*N)

          IWORK is INTEGER array, dimension (N)

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

 STBTRS solves a triangular system of the form

    A * X = B  or  A**T * X = B,

 where A is a triangular band matrix of order N, and B is an
 N-by NRHS matrix.  A check is made to verify that A is nonsingular.

          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.

          TRANS is CHARACTER*1
          Specifies the form the system of equations:
          = 'N':  A * X = B  (No transpose)
          = 'T':  A**T * X = B  (Transpose)
          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)

          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.

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

          KD is INTEGER
          The number of superdiagonals or subdiagonals of the
          triangular band matrix A.  KD >= 0.

          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrix B.  NRHS >= 0.

          AB is REAL array, dimension (LDAB,N)
          The upper or lower triangular band matrix A, stored in the
          first kd+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(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).
          If DIAG = 'U', the diagonal elements of A are not referenced
          and are assumed to be 1.

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

          B is REAL array, dimension (LDB,NRHS)
          On entry, the right hand side matrix B.
          On exit, if INFO = 0, the solution matrix X.

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

          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the i-th diagonal element of A is zero,
                indicating that the matrix is singular and the
                solutions X have not been computed.

 Level 3 BLAS like routine for A in RFP Format.

 STFSM  solves the matrix equation

    op( A )*X = alpha*B  or  X*op( A ) = alpha*B

 where alpha is a scalar, X and B are m by n matrices, A is a unit, or
 non-unit,  upper or lower triangular matrix  and  op( A )  is one  of

    op( A ) = A   or   op( A ) = A**T.

 A is in Rectangular Full Packed (RFP) Format.

 The matrix X is overwritten on B.

          TRANSR is CHARACTER*1
          = 'N':  The Normal Form of RFP A is stored;
          = 'T':  The Transpose Form of RFP A is stored.

          SIDE is CHARACTER*1
           On entry, SIDE specifies whether op( A ) appears on the left
           or right of X as follows:

              SIDE = 'L' or 'l'   op( A )*X = alpha*B.

              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.

           Unchanged on exit.

          UPLO is CHARACTER*1
           On entry, UPLO specifies whether the RFP matrix A came from
           an upper or lower triangular matrix as follows:
           UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
           UPLO = 'L' or 'l' RFP A came from a  lower triangular matrix

           Unchanged on exit.

          TRANS is CHARACTER*1
           On entry, TRANS  specifies the form of op( A ) to be used
           in the matrix multiplication as follows:

              TRANS  = 'N' or 'n'   op( A ) = A.

              TRANS  = 'T' or 't'   op( A ) = A'.

           Unchanged on exit.

          DIAG is CHARACTER*1
           On entry, DIAG specifies whether or not RFP A is unit
           triangular as follows:

              DIAG = 'U' or 'u'   A is assumed to be unit triangular.

              DIAG = 'N' or 'n'   A is not assumed to be unit
                                  triangular.

           Unchanged on exit.

          M is INTEGER
           On entry, M specifies the number of rows of B. M must be at
           least zero.
           Unchanged on exit.

          N is INTEGER
           On entry, N specifies the number of columns of B.  N must be
           at least zero.
           Unchanged on exit.

          ALPHA is REAL
           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
           zero then  A is not referenced and  B need not be set before
           entry.
           Unchanged on exit.

          A is REAL array, dimension (NT)
           NT = N*(N+1)/2. On entry, the matrix A in RFP Format.
           RFP Format is described by TRANSR, UPLO and N as follows:
           If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
           K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If
           TRANSR = 'T' then RFP is the transpose of RFP A as
           defined when TRANSR = 'N'. The contents of RFP A are defined
           by UPLO as follows: If UPLO = 'U' the RFP A contains the NT
           elements of upper packed A either in normal or
           transpose Format. If UPLO = 'L' the RFP A contains
           the NT elements of lower packed A either in normal or
           transpose Format. The LDA of RFP A is (N+1)/2 when
           TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
           even and is N when is odd.
           See the Note below for more details. Unchanged on exit.

          B is REAL array, DIMENSION (LDB,N)
           Before entry,  the leading  m by n part of the array  B must
           contain  the  right-hand  side  matrix  B,  and  on exit  is
           overwritten by the solution matrix  X.

          LDB is INTEGER
           On entry, LDB specifies the first dimension of B as declared
           in  the  calling  (sub)  program.   LDB  must  be  at  least
           max( 1, m ).
           Unchanged on exit.

  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

 STFTRI computes the inverse of a triangular matrix A stored in RFP
 format.

 This is a Level 3 BLAS version of the algorithm.

          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.

          UPLO is CHARACTER*1
          = 'U':  A is upper triangular;
          = 'L':  A is lower triangular.

          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.

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

          A is REAL array, dimension (NT);
          NT=N*(N+1)/2. On entry, the triangular factor of a Hermitian
          Positive Definite matrix A in RFP format. RFP format is
          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
          the transpose of RFP A as defined when
          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
          follows: If UPLO = 'U' the RFP A contains the nt elements of
          upper packed A; If UPLO = 'L' the RFP A contains the nt
          elements of lower packed A. The LDA of RFP A is (N+1)/2 when
          TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
          even and N is odd. See the Note below for more details.

          On exit, the (triangular) inverse of the original matrix, in
          the same storage format.

          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
               matrix is singular and its inverse can not be computed.