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

          JOBU1 is CHARACTER


          = 'Y':      U1 is updated;


          otherwise:  U1 is not updated.

JOBU2

          JOBU2 is CHARACTER


          = 'Y':      U2 is updated;


          otherwise:  U2 is not updated.

JOBV1T

          JOBV1T is CHARACTER


          = 'Y':      V1T is updated;


          otherwise:  V1T is not updated.

JOBV2T

          JOBV2T is CHARACTER


          = 'Y':      V2T is updated;


          otherwise:  V2T is not updated.

TRANS

          TRANS is CHARACTER


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


                      order;


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


                      major order.

M

          M is INTEGER


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


          bidiagonal-block form.

P

          P is INTEGER


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

Q

          Q is INTEGER


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


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

THETA

          THETA is 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

          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

          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

          LDU1 is INTEGER


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

U2

          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

          LDU2 is INTEGER


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

V1T

          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

          LDV1T is INTEGER


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

V2T

          V2T is REAL array, dimension (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

          LDV2T is INTEGER


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

B11D

          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

          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

          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

          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

          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

          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

          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

          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

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= MAX(1,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

          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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          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

          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

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


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


          reduced.  It is assumed that A is already upper triangular


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


          normally set by a previous call to 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

          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

          LDA is INTEGER


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

B

          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

          LDB is INTEGER


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

Q

          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

          LDQ is INTEGER


          The leading dimension of the array Q.


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

Z

          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

          LDZ is INTEGER


          The leading dimension of the array Z.


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

WORK

          WORK is REAL array, dimension (LWORK)


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

LWORK

          LWORK is INTEGER


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


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


          optimal blocksize.


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


          only calculates the optimal size of the WORK array, returns


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


          message related to LWORK is issued by XERBLA.

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


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


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


  (BIT 2008).

 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

          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

          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

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


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


          reduced.  It is assumed that A is already upper triangular


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


          normally set by a previous call to 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

          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

          LDA is INTEGER


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

B

          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

          LDB is INTEGER


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

Q

          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

          LDQ is INTEGER


          The leading dimension of the array Q.


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

Z

          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

          LDZ is INTEGER


          The leading dimension of the array Z.


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

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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


  an unblocked reduction, as described in _Matrix_Computations_,


  by Golub and Van Loan (Johns Hopkins Press.)

 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

          N is INTEGER


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

M

          M is INTEGER


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

P

          P is INTEGER


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

A

          A is 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

          LDA is INTEGER


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

TAUA

          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

          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

          LDB is INTEGER


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

TAUB

          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

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


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


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


          where NB1 is the optimal blocksize for the QR factorization


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


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


          blocksize for a call of 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The matrix Q is represented as a product of elementary reflectors


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


  Each H(i) has the form


     H(i) = I - taua * v * v**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

          M is INTEGER


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

P

          P is INTEGER


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

N

          N is INTEGER


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

A

          A is 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

          LDA is INTEGER


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

TAUA

          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

          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

          LDB is INTEGER


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

TAUB

          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

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


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


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


          where NB1 is the optimal blocksize for the RQ factorization


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


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


          blocksize for a call of 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The matrix Q is represented as a product of elementary reflectors


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


  Each H(i) has the form


     H(i) = I - taua * v * v**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

          JOBU is CHARACTER*1


          = 'U':  Orthogonal matrix U is computed;


          = 'N':  U is not computed.

JOBV

          JOBV is CHARACTER*1


          = 'V':  Orthogonal matrix V is computed;


          = 'N':  V is not computed.

JOBQ

          JOBQ is CHARACTER*1


          = 'Q':  Orthogonal matrix Q is computed;


          = 'N':  Q is not computed.

M

          M is INTEGER


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

P

          P is INTEGER


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

N

          N is INTEGER


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

A

          A is 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

          LDA is INTEGER


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

B

          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

          LDB is INTEGER


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

TOLA

          TOLA is REAL

TOLB

          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

          K is INTEGER

L

          L is INTEGER


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


          described in Purpose section.


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

U

          U is REAL array, dimension (LDU,M)


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


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

LDU

          LDU is INTEGER


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


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

V

          V is REAL array, dimension (LDV,P)


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


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

LDV

          LDV is INTEGER


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


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

Q

          Q is REAL array, dimension (LDQ,N)


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


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

LDQ

          LDQ is INTEGER


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


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

IWORK

          IWORK is INTEGER array, dimension (N)

TAU

          TAU is REAL array, dimension (N)

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


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


          only calculates the optimal size of the WORK array, returns


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


          message related to LWORK is issued by XERBLA.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  The subroutine uses LAPACK subroutine 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

          JOBV is CHARACTER*1


          Specifies whether the output from this procedure is used


          to compute the matrix V:


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


                 by postmulyiplying the N-by-N array V.


                (See the description of V.)


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


                 by postmulyiplying the MV-by-N array V.


                (See the descriptions of MV and V.)


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

M

          M is INTEGER


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

N

          N is INTEGER


          The number of columns of the input matrix A.


          M >= N >= 0.

A

          A is 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

          LDA is INTEGER


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

D

          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

          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

          MV is INTEGER


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


                           sequence of Jacobi rotations.


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

V

          V is REAL array, dimension (LDV,N)


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


                           sequence of Jacobi rotations.


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


                           sequence of Jacobi rotations.


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

LDV

          LDV is INTEGER


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


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


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

EPS

          EPS is REAL


          EPS = SLAMCH('Epsilon')

SFMIN

          SFMIN is REAL


          SFMIN = SLAMCH('Safe Minimum')

TOL

          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)))) > TOL.

NSWEEP

          NSWEEP is INTEGER


          NSWEEP is the number of sweeps of Jacobi rotations to be


          performed.

WORK

          WORK is REAL array, dimension (LWORK)

LWORK

          LWORK is INTEGER


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

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

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

 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 tuning 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 tuning parameter.


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


 is given in TOL.

JOBV

          JOBV is CHARACTER*1


          Specifies whether the output from this procedure is used


          to compute the matrix V:


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


                 by postmulyiplying the N-by-N array V.


                (See the description of V.)


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


                 by postmulyiplying the MV-by-N array V.


                (See the descriptions of MV and V.)


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

M

          M is INTEGER


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

N

          N is INTEGER


          The number of columns of the input matrix A.


          M >= N >= 0.

N1

          N1 is INTEGER


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


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

A

          A is 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

          LDA is INTEGER


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

D

          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

          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

          MV is INTEGER


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


                           sequence of Jacobi rotations.


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

V

          V is REAL array, dimension (LDV,N)


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


                           sequence of Jacobi rotations.


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


                           sequence of Jacobi rotations.


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

LDV

          LDV is INTEGER


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


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


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

EPS

          EPS is REAL


          EPS = SLAMCH('Epsilon')

SFMIN

          SFMIN is REAL


          SFMIN = SLAMCH('Safe Minimum')

TOL

          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)))) > TOL.

NSWEEP

          NSWEEP is INTEGER


          NSWEEP is the number of sweeps of Jacobi rotations to be


          performed.

WORK

          WORK is REAL array, dimension (LWORK)

LWORK

          LWORK is INTEGER


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

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

 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

          SIDE is CHARACTER*1


          = 'R': compute right eigenvectors only;


          = 'L': compute left eigenvectors only;


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

EIGSRC

          EIGSRC is CHARACTER*1


          Specifies the source of eigenvalues supplied in (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

          INITV is CHARACTER*1


          = 'N': no initial vectors are supplied;


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


                 VL and/or VR.

SELECT

          SELECT is LOGICAL array, dimension (N)


          Specifies the eigenvectors to be computed. To select the


          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

          N is INTEGER


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

H

          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

          LDH is INTEGER


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

WR

          WR is REAL array, dimension (N)

WI

          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

          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

          LDVL is INTEGER


          The leading dimension of the array VL.


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

VR

          VR is 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

          LDVR is INTEGER


          The leading dimension of the array VR.


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

MM

          MM is INTEGER


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

M

          M is INTEGER


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


          store the eigenvectors; each selected real eigenvector


          occupies one column and each selected complex eigenvector


          occupies two columns.

WORK

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

IFAILL

          IFAILL is INTEGER array, dimension (MM)


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


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


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


          eigenvector converged satisfactorily. If 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

          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

          INFO is INTEGER


          = 0:  successful exit


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


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


                failed to converge; see IFAILL and IFAILR for further


                details.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  Each eigenvector is normalized so that the element of largest


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


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

    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

          JOB is CHARACTER*1


           = 'E':  compute eigenvalues only;


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

COMPZ

          COMPZ is CHARACTER*1


           = 'N':  no Schur vectors are computed;


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


                   of Schur vectors of H is returned;


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


                   the product Q*Z is returned.

N

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          IHI is INTEGER


           It is assumed that H is already upper triangular in rows


           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally


           set by a previous call to 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 > 0, then 1 <= ILO <= IHI <= N.


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

H

          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) < 0. If INFO = 0 and JOB = 'E', the


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


           H when INFO > 0 is given under the description of INFO


           below.)


           Unlike earlier versions of SHSEQR, this subroutine may


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


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

LDH

          LDH is INTEGER


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

WR

          WR is REAL array, dimension (N)

WI

          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) > 0 and


           WI(i+1) < 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

          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 > 0 is given under


           the description of INFO below.)

LDZ

          LDZ is INTEGER


           The leading dimension of the array Z.  if COMPZ = 'I' or


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

WORK

          WORK is REAL array, dimension (LWORK)


           On exit, if INFO = 0, WORK(1) returns an estimate of


           the optimal value for LWORK.

LWORK

          LWORK is INTEGER


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


           is sufficient and delivers very good and sometimes


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


           may be required for optimal performance.  A workspace


           query is recommended to determine the optimal workspace


           size.


           If LWORK = -1, then 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

          INFO is INTEGER


             = 0:  successful exit


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


                    value


             > 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 > 0 and JOB = 'E', then on exit, the


                remaining unconverged eigenvalues are the eigen-


                values of the upper Hessenberg matrix rows and


                columns ILO through INFO of the final, output


                value of H.


                If INFO > 0 and JOB   = 'S', then on exit


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


                where U is an orthogonal matrix.  The final


                value of H is upper Hessenberg and quasi-triangular


                in rows and columns INFO+1 through IHI.


                If INFO > 0 and 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 > 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 > 0 and COMPZ = 'N', then Z is not


                accessed.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

             Default values supplied by


             ILAENV(ISPEC,'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) <= 500 is NS.


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


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


                      details.)  Default: 14% of deflation window


                      size.


            ISPEC=15: Number of simultaneous shifts in a multishift


                      QR iteration.


                      If IHI-ILO+1 is ...


                      greater than      ...but less    ... the


                      or equal to ...      than        default is


                           1               30          NS =   2(+)


                          30               60          NS =   4(+)


                          60              150          NS =  10(+)


                         150              590          NS =  **


                         590             3000          NS =  64


                        3000             6000          NS = 128


                        6000             infinity      NS = 256


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


                       are passed to the implicit double shift routine


                       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.

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


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


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


  929--947, 2002.

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

    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

          N is INTEGER


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


     matrix A.  N >= 0.

NZ

          NZ is INTEGER


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


     guard against spuriously zero residuals. Default value is N.

NRHS

          NRHS is INTEGER


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


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

RES

          RES is REAL array, dimension (N,NRHS)


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


     error formula above.

AYB

          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

          BERR is REAL array, dimension (NRHS)


     The componentwise relative backward error from the formula above.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

    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


    arithmetic, but not for decimal.

N

          N is INTEGER


            The length of vectors X, Y, and W.

X

          X is REAL array, dimension (N)


            The first part of the doubled-single accumulation vector.

Y

          Y is REAL array, dimension (N)


            The second part of the doubled-single accumulation vector.

W

          W is REAL array, dimension (N)


            The vector to be added.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          ICOMPQ is INTEGER


         Specifies whether singular vectors are to be computed in


         factored form:


         = 0: Left singular vector matrix.


         = 1: Right singular vector matrix.

NL

          NL is INTEGER


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

NR

          NR is INTEGER


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

SQRE

          SQRE is INTEGER


         = 0: the lower block is an NR-by-NR square matrix.


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


         The bidiagonal matrix has row dimension N = NL + NR + 1,


         and column dimension M = N + SQRE.

NRHS

          NRHS is INTEGER


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

B

          B is 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

          LDB is INTEGER


         The leading dimension of B. LDB must be at least


         max(1,MAX( M, N ) ).

BX

          BX is REAL array, dimension ( LDBX, NRHS )

LDBX

          LDBX is INTEGER


         The leading dimension of BX.

PERM

          PERM is INTEGER array, dimension ( N )


         The permutations (from deflation and sorting) applied


         to the two blocks.

GIVPTR

          GIVPTR is INTEGER


         The number of Givens rotations which took place in this


         subproblem.

GIVCOL

          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )


         Each pair of numbers indicates a pair of rows/columns


         involved in a Givens rotation.

LDGCOL

          LDGCOL is INTEGER


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

GIVNUM

          GIVNUM is REAL array, dimension ( LDGNUM, 2 )


         Each number indicates the C or S value used in the


         corresponding Givens rotation.

LDGNUM

          LDGNUM is INTEGER


         The leading dimension of arrays DIFR, POLES and


         GIVNUM, must be at least K.

POLES

          POLES is 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

          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

          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

          Z is REAL array, dimension ( K )


         Contain the components of the deflation-adjusted updating row


         vector.

K

          K is INTEGER


         Contains the dimension of the non-deflated matrix,


         This is the order of the related secular equation. 1 <= K <=N.

C

          C is 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

          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

          WORK is REAL array, dimension ( K )

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

 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

          ICOMPQ is INTEGER


         Specifies whether the left or the right singular vector


         matrix is involved.


         = 0: Left singular vector matrix


         = 1: Right singular vector matrix

SMLSIZ

          SMLSIZ is INTEGER


         The maximum size of the subproblems at the bottom of the


         computation tree.

N

          N is INTEGER


         The row and column dimensions of the upper bidiagonal matrix.

NRHS

          NRHS is INTEGER


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

B

          B is 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

          LDB is INTEGER


         The leading dimension of B in the calling subprogram.


         LDB must be at least max(1,MAX( M, N ) ).

BX

          BX is REAL array, dimension ( LDBX, NRHS )


         On exit, the result of applying the left or right singular


         vector matrix to B.

LDBX

          LDBX is INTEGER


         The leading dimension of BX.

U

          U is REAL array, dimension ( LDU, SMLSIZ ).


         On entry, U contains the left singular vector matrices of all


         subproblems at the bottom level.

LDU

          LDU is INTEGER, LDU = > N.


         The leading dimension of arrays U, VT, DIFL, DIFR,


         POLES, GIVNUM, and Z.

VT

          VT is REAL array, dimension ( LDU, SMLSIZ+1 ).


         On entry, VT**T contains the right singular vector matrices of


         all subproblems at the bottom level.

K

          K is INTEGER array, dimension ( N ).

DIFL

          DIFL is REAL array, dimension ( LDU, NLVL ).


         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.

DIFR

          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

          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

          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

          GIVPTR is INTEGER array, dimension ( N ).


         On entry, GIVPTR( I ) records the number of Givens


         rotations performed on the I-th problem on the computation


         tree.

GIVCOL

          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).


         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the


         locations of Givens rotations performed on the I-th level on


         the computation tree.

LDGCOL

          LDGCOL is INTEGER, LDGCOL = > N.


         The leading dimension of arrays GIVCOL and PERM.

PERM

          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).


         On entry, PERM(*, I) records permutations done on the I-th


         level of the computation tree.

GIVNUM

          GIVNUM is 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

          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

          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

          WORK is REAL array, dimension (N)

IWORK

          IWORK is INTEGER array, dimension (3*N)

INFO

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

 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

          UPLO is CHARACTER*1


         = 'U': D and E define an upper bidiagonal matrix.


         = 'L': D and E define a  lower bidiagonal matrix.

SMLSIZ

          SMLSIZ is INTEGER


         The maximum size of the subproblems at the bottom of the


         computation tree.

N

          N is INTEGER


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

NRHS

          NRHS is INTEGER


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

D

          D is 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

          E is REAL array, dimension (N-1)


         Contains the super-diagonal entries of the bidiagonal matrix.


         On exit, E has been destroyed.

B

          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

          LDB is INTEGER


         The leading dimension of B in the calling subprogram.


         LDB must be at least max(1,N).

RCOND

          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

          RANK is INTEGER


         The number of singular values of A greater than RCOND times


         the largest singular value.

WORK

          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

          IWORK is INTEGER array, dimension at least


         (3*N*NLVL + 11*N)

INFO

          INFO is INTEGER


         = 0:  successful exit.


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


         > 0:  The algorithm failed to compute a singular value while


               working on the submatrix lying in rows and columns


               INFO/(N+1) through MOD(INFO,N+1).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

 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

    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

          NORM is CHARACTER*1


          Specifies the value to be returned in SLANSF as described


          above.

TRANSR

          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

          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

          N is INTEGER


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


          set to zero.

A

          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

          WORK is REAL array, dimension (MAX(1,LWORK)),


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


          WORK is not referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  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 a matrix:


   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

          M is INTEGER


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

N

          N is INTEGER


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

D

          D is REAL array, length M


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

X

          X is REAL array, dimension (LDX,N)


     On entry, the matrix X to be scaled by D.


     On exit, the scaled matrix.

LDX

          LDX is INTEGER


     The leading dimension of the matrix X. LDX >= M.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': form  H * C


          = 'R': form  C * H

M

          M is INTEGER


          The number of rows of the matrix C.

N

          N is INTEGER


          The number of columns of the matrix C.

L

          L is INTEGER


          The number of entries of the vector V containing


          the meaningful part of the Householder vectors.


          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

V

          V is 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

          INCV is INTEGER


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

TAU

          TAU is REAL


          The value tau in the representation of H.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                         (N) if SIDE = 'L'


                      or (M) if SIDE = 'R'

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

 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

          SIDE is CHARACTER*1


          = 'L': apply H or H**T from the Left


          = 'R': apply H or H**T from the Right

TRANS

          TRANS is CHARACTER*1


          = 'N': apply H (No transpose)


          = 'C': apply H**T (Transpose)

DIRECT

          DIRECT is CHARACTER*1


          Indicates how H is formed from a product of elementary


          reflectors


          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)


          = 'B': H = H(k) . . . H(2) H(1) (Backward)

STOREV

          STOREV is CHARACTER*1


          Indicates how the vectors which define the elementary


          reflectors are stored:


          = 'C': Columnwise                        (not supported yet)


          = 'R': Rowwise

M

          M is INTEGER


          The number of rows of the matrix C.

N

          N is INTEGER


          The number of columns of the matrix C.

K

          K is INTEGER


          The order of the matrix T (= the number of elementary


          reflectors whose product defines the block reflector).

L

          L is INTEGER


          The number of columns of the matrix V containing the


          meaningful part of the Householder reflectors.


          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

V

          V is REAL array, dimension (LDV,NV).


          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.

T

          T is REAL array, dimension (LDT,K)


          The triangular K-by-K matrix T in the representation of the


          block reflector.

LDT

          LDT is INTEGER


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

C

          C is 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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (LDWORK,K)

LDWORK

          LDWORK is INTEGER


          The leading dimension of the array WORK.


          If SIDE = 'L', LDWORK >= max(1,N);


          if SIDE = 'R', LDWORK >= max(1,M).

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

 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

          DIRECT is CHARACTER*1


          Specifies the order in which the elementary reflectors are


          multiplied to form the block reflector:


          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)


          = 'B': H = H(k) . . . H(2) H(1) (Backward)

STOREV

          STOREV is CHARACTER*1


          Specifies how the vectors which define the elementary


          reflectors are stored (see also Further Details):


          = 'C': columnwise                        (not supported yet)


          = 'R': rowwise

N

          N is INTEGER


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

K

          K is INTEGER


          The order of the triangular factor T (= the number of


          elementary reflectors). K >= 1.

V

          V is REAL array, dimension


                               (LDV,K) if STOREV = 'C'


                               (LDV,N) if STOREV = 'R'


          The matrix V. See further details.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i).

T

          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

          LDT is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

  The shape of the matrix V and the storage of the vectors which define


  the H(i) is best illustrated by the following example with n = 5 and


  k = 3. The elements equal to 1 are not stored; the corresponding


  array elements are modified but restored on exit. The rest of the


  array is not used.


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


                                              ______V_____


         ( v1 v2 v3 )                        /            \


         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )


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


         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )


         ( v1 v2 v3 )


            .  .  .


            .  .  .


            1  .  .


               1  .


                  1


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


                                                        ______V_____


            1                                          /            \


            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )


            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )


            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )


            .  .  .


         ( v1 v2 v3 )


         ( v1 v2 v3 )


     V = ( v1 v2 v3 )


         ( v1 v2 v3 )


         ( v1 v2 v3 )

 SLASCL2 performs a diagonal scaling on a matrix:


   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

          M is INTEGER


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

N

          N is INTEGER


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

D

          D is REAL array, length M


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

X

          X is REAL array, dimension (LDX,N)


     On entry, the matrix X to be scaled by D.


     On exit, the scaled matrix.

LDX

          LDX is INTEGER


     The leading dimension of the matrix X. LDX >= M.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

L

          L is INTEGER


          The number of columns of the matrix A containing the


          meaningful part of the Householder vectors. N-M >= L >= 0.

A

          A is 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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (M)


          The scalar factors of the elementary reflectors.

WORK

          WORK is REAL array, dimension (M)

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

  The factorization is obtained by Householder's method.  The kth


  transformation matrix, Z( k ), which is used to introduce zeros into


  the ( m - k + 1 )th row of A, is given in the form


     Z( k ) = ( I     0   ),


              ( 0  T( k ) )


  where


     T( k ) = I - tau*u( k )*u( k )**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

          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

          N is INTEGER


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

AP

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


          The vectors which define the elementary reflectors, as


          returned by SSPTRD.

TAU

          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

          Q is REAL array, dimension (LDQ,N)


          The N-by-N orthogonal matrix Q.

LDQ

          LDQ is INTEGER


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

WORK

          WORK is REAL array, dimension (N-1)

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

UPLO

          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

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

AP

          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

          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

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                                   (N) if SIDE = 'L'


                                   (M) if SIDE = 'R'

INFO

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          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

          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

          M is INTEGER


          The number of rows and columns in X.

P

          P is INTEGER


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

Q

          Q is INTEGER


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


          MIN(P,M-P,M-Q).

X11

          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

          LDX11 is INTEGER


          The leading dimension of X11. If TRANS = 'N', then LDX11 >=


          P; else LDX11 >= Q.

X12

          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

          LDX12 is INTEGER


          The leading dimension of X12. If TRANS = 'N', then LDX12 >=


          P; else LDX11 >= M-Q.

X21

          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

          LDX21 is INTEGER


          The leading dimension of X21. If TRANS = 'N', then LDX21 >=


          M-P; else LDX21 >= Q.

X22

          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

          LDX22 is INTEGER


          The leading dimension of X22. If TRANS = 'N', then LDX22 >=


          M-P; else LDX22 >= M-Q.

THETA

          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

          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

          TAUP1 is REAL array, dimension (P)


          The scalar factors of the elementary reflectors that define


          P1.

TAUP2

          TAUP2 is REAL array, dimension (M-P)


          The scalar factors of the elementary reflectors that define


          P2.

TAUQ1

          TAUQ1 is REAL array, dimension (Q)


          The scalar factors of the elementary reflectors that define


          Q1.

TAUQ2

          TAUQ2 is REAL array, dimension (M-Q)


          The scalar factors of the elementary reflectors that define


          Q2.

WORK

          WORK is REAL array, dimension (LWORK)

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

          INFO is INTEGER


          = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

 SORBDB1 simultaneously bidiagonalizes the blocks of a tall and skinny


 matrix X with orthonomal columns:


                            [ B11 ]


      [ X11 ]   [ P1 |    ] [  0  ]


      [-----] = [---------] [-----] Q1**T .


      [ X21 ]   [    | P2 ] [ B21 ]


                            [  0  ]


 X11 is P-by-Q, and X21 is (M-P)-by-Q. Q must be no larger than P,


 M-P, or M-Q. Routines SORBDB2, SORBDB3, and SORBDB4 handle cases in


 which Q is not the minimum dimension.


 The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),


 and (M-Q)-by-(M-Q), respectively. They are represented implicitly by


 Householder vectors.


 B11 and B12 are Q-by-Q bidiagonal matrices represented implicitly by


 angles THETA, PHI.

M

          M is INTEGER


           The number of rows X11 plus the number of rows in X21.

P

          P is INTEGER


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

Q

          Q is INTEGER


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


           MIN(P,M-P,M-Q).

X11

          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

          LDX11 is INTEGER


           The leading dimension of X11. LDX11 >= P.

X21

          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

          LDX21 is INTEGER


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

THETA

          THETA is REAL array, dimension (Q)


           The entries of the bidiagonal blocks B11, B21 are defined by


           THETA and PHI. See Further Details.

PHI

          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

          TAUP1 is REAL array, dimension (P)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is REAL array, dimension (M-P)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is REAL array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

WORK

          WORK is REAL array, dimension (LWORK)

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

          INFO is INTEGER


           = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

 SORBDB2 simultaneously bidiagonalizes the blocks of a tall and skinny


 matrix X with orthonomal columns:


                            [ B11 ]


      [ X11 ]   [ P1 |    ] [  0  ]


      [-----] = [---------] [-----] Q1**T .


      [ X21 ]   [    | P2 ] [ B21 ]


                            [  0  ]


 X11 is P-by-Q, and X21 is (M-P)-by-Q. P must be no larger than M-P,


 Q, or M-Q. Routines SORBDB1, SORBDB3, and SORBDB4 handle cases in


 which P is not the minimum dimension.


 The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),


 and (M-Q)-by-(M-Q), respectively. They are represented implicitly by


 Householder vectors.


 B11 and B12 are P-by-P bidiagonal matrices represented implicitly by


 angles THETA, PHI.

M

          M is INTEGER


           The number of rows X11 plus the number of rows in X21.

P

          P is INTEGER


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

Q

          Q is INTEGER


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

X11

          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

          LDX11 is INTEGER


           The leading dimension of X11. LDX11 >= P.

X21

          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

          LDX21 is INTEGER


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

THETA

          THETA is REAL array, dimension (Q)


           The entries of the bidiagonal blocks B11, B21 are defined by


           THETA and PHI. See Further Details.

PHI

          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

          TAUP1 is REAL array, dimension (P-1)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is REAL array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is REAL array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

WORK

          WORK is REAL array, dimension (LWORK)

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

          INFO is INTEGER


           = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

 SORBDB3 simultaneously bidiagonalizes the blocks of a tall and skinny


 matrix X with orthonomal columns:


                            [ B11 ]


      [ X11 ]   [ P1 |    ] [  0  ]


      [-----] = [---------] [-----] Q1**T .


      [ X21 ]   [    | P2 ] [ B21 ]


                            [  0  ]


 X11 is P-by-Q, and X21 is (M-P)-by-Q. M-P must be no larger than P,


 Q, or M-Q. Routines SORBDB1, SORBDB2, and SORBDB4 handle cases in


 which M-P is not the minimum dimension.


 The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),


 and (M-Q)-by-(M-Q), respectively. They are represented implicitly by


 Householder vectors.


 B11 and B12 are (M-P)-by-(M-P) bidiagonal matrices represented


 implicitly by angles THETA, PHI.

M

          M is INTEGER


           The number of rows X11 plus the number of rows in X21.

P

          P is INTEGER


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

Q

          Q is INTEGER


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

X11

          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

          LDX11 is INTEGER


           The leading dimension of X11. LDX11 >= P.

X21

          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

          LDX21 is INTEGER


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

THETA

          THETA is REAL array, dimension (Q)


           The entries of the bidiagonal blocks B11, B21 are defined by


           THETA and PHI. See Further Details.

PHI

          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

          TAUP1 is REAL array, dimension (P)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is REAL array, dimension (M-P)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is REAL array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

WORK

          WORK is REAL array, dimension (LWORK)

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

          INFO is INTEGER


           = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

 SORBDB4 simultaneously bidiagonalizes the blocks of a tall and skinny


 matrix X with orthonomal columns:


                            [ B11 ]


      [ X11 ]   [ P1 |    ] [  0  ]


      [-----] = [---------] [-----] Q1**T .


      [ X21 ]   [    | P2 ] [ B21 ]


                            [  0  ]


 X11 is P-by-Q, and X21 is (M-P)-by-Q. M-Q must be no larger than P,


 M-P, or Q. Routines SORBDB1, SORBDB2, and SORBDB3 handle cases in


 which M-Q is not the minimum dimension.


 The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),


 and (M-Q)-by-(M-Q), respectively. They are represented implicitly by


 Householder vectors.


 B11 and B12 are (M-Q)-by-(M-Q) bidiagonal matrices represented


 implicitly by angles THETA, PHI.

M

          M is INTEGER


           The number of rows X11 plus the number of rows in X21.

P

          P is INTEGER


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

Q

          Q is INTEGER


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


           M-Q <= min(P,M-P,Q).

X11

          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

          LDX11 is INTEGER


           The leading dimension of X11. LDX11 >= P.

X21

          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

          LDX21 is INTEGER


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

THETA

          THETA is REAL array, dimension (Q)


           The entries of the bidiagonal blocks B11, B21 are defined by


           THETA and PHI. See Further Details.

PHI

          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

          TAUP1 is REAL array, dimension (M-Q)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is REAL array, dimension (M-Q)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is REAL array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

PHANTOM

          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

          WORK is REAL array, dimension (LWORK)

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

          INFO is INTEGER


           = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

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

 SORBDB5 orthogonalizes the column vector


      X = [ X1 ]


          [ X2 ]


 with respect to the columns of


      Q = [ Q1 ] .


          [ Q2 ]


 The columns of Q must be orthonormal.


 If the projection is zero according to Kahan's 'twice is enough'


 criterion, then some other vector from the orthogonal complement


 is returned. This vector is chosen in an arbitrary but deterministic


 way.

M1

          M1 is INTEGER


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

M2

          M2 is INTEGER


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

N

          N is INTEGER


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

X1

          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

          INCX1 is INTEGER


           Increment for entries of X1.

X2

          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

          INCX2 is INTEGER


           Increment for entries of X2.

Q1

          Q1 is REAL array, dimension (LDQ1, N)


           The top part of the orthonormal basis matrix.

LDQ1

          LDQ1 is INTEGER


           The leading dimension of Q1. LDQ1 >= M1.

Q2

          Q2 is REAL array, dimension (LDQ2, N)


           The bottom part of the orthonormal basis matrix.

LDQ2

          LDQ2 is INTEGER


           The leading dimension of Q2. LDQ2 >= M2.

WORK

          WORK is REAL array, dimension (LWORK)

LWORK

          LWORK is INTEGER


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

INFO

          INFO is INTEGER


           = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 SORBDB6 orthogonalizes the column vector


      X = [ X1 ]


          [ X2 ]


 with respect to the columns of


      Q = [ Q1 ] .


          [ Q2 ]


 The Euclidean norm of X must be one and the columns of Q must be


 orthonormal. The orthogonalized vector will be zero if and only if it


 lies entirely in the range of Q.


 The projection is computed with at most two iterations of the


 classical Gram-Schmidt algorithm, see


 * L. Giraud, J. Langou, M. Rozložník. 'On the round-off error


   analysis of the Gram-Schmidt algorithm with reorthogonalization.'


   2002. CERFACS Technical Report No. TR/PA/02/33. URL:


   https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf

M1

          M1 is INTEGER


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

M2

          M2 is INTEGER


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

N

          N is INTEGER


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

X1

          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

          INCX1 is INTEGER


           Increment for entries of X1.

X2

          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

          INCX2 is INTEGER


           Increment for entries of X2.

Q1

          Q1 is REAL array, dimension (LDQ1, N)


           The top part of the orthonormal basis matrix.

LDQ1

          LDQ1 is INTEGER


           The leading dimension of Q1. LDQ1 >= M1.

Q2

          Q2 is REAL array, dimension (LDQ2, N)


           The bottom part of the orthonormal basis matrix.

LDQ2

          LDQ2 is INTEGER


           The leading dimension of Q2. LDQ2 >= M2.

WORK

          WORK is REAL array, dimension (LWORK)

LWORK

          LWORK is INTEGER


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

INFO

          INFO is INTEGER


           = 0:  successful exit.


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          JOBU1 is CHARACTER


          = 'Y':      U1 is computed;


          otherwise:  U1 is not computed.

JOBU2

          JOBU2 is CHARACTER


          = 'Y':      U2 is computed;


          otherwise:  U2 is not computed.

JOBV1T

          JOBV1T is CHARACTER


          = 'Y':      V1T is computed;


          otherwise:  V1T is not computed.

JOBV2T

          JOBV2T is CHARACTER


          = 'Y':      V2T is computed;


          otherwise:  V2T is not computed.

TRANS

          TRANS is CHARACTER


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


                      order;


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


                      major order.

SIGNS

          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

          M is INTEGER


          The number of rows and columns in X.

P

          P is INTEGER


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

Q

          Q is INTEGER


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

X11

          X11 is REAL array, dimension (LDX11,Q)


          On entry, part of the orthogonal matrix whose CSD is desired.

LDX11

          LDX11 is INTEGER


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

X12

          X12 is REAL array, dimension (LDX12,M-Q)


          On entry, part of the orthogonal matrix whose CSD is desired.

LDX12

          LDX12 is INTEGER


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

X21

          X21 is REAL array, dimension (LDX21,Q)


          On entry, part of the orthogonal matrix whose CSD is desired.

LDX21

          LDX21 is INTEGER


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

X22

          X22 is REAL array, dimension (LDX22,M-Q)


          On entry, part of the orthogonal matrix whose CSD is desired.

LDX22

          LDX22 is INTEGER


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

THETA

          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

          U1 is REAL array, dimension (LDU1,P)


          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.

LDU1

          LDU1 is INTEGER


          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=


          MAX(1,P).

U2

          U2 is REAL array, dimension (LDU2,M-P)


          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal


          matrix U2.

LDU2

          LDU2 is INTEGER


          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=


          MAX(1,M-P).

V1T

          V1T is REAL array, dimension (LDV1T,Q)


          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal


          matrix V1**T.

LDV1T

          LDV1T is INTEGER


          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=


          MAX(1,Q).

V2T

          V2T is REAL array, dimension (LDV2T,M-Q)


          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal


          matrix V2**T.

LDV2T

          LDV2T is INTEGER


          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=


          MAX(1,M-Q).

WORK

          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

          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

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

INFO

          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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 SORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with


 orthonormal columns that has been partitioned into a 2-by-1 block


 structure:


                                [  I1 0  0 ]


                                [  0  C  0 ]


          [ X11 ]   [ U1 |    ] [  0  0  0 ]


      X = [-----] = [---------] [----------] V1**T .


          [ X21 ]   [    | U2 ] [  0  0  0 ]


                                [  0  S  0 ]


                                [  0  0  I2]


 X11 is P-by-Q. The orthogonal matrices U1, U2, and V1 are P-by-P,


 (M-P)-by-(M-P), and Q-by-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). I1 is a K1-by-K1 identity matrix and I2 is a


 K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).

JOBU1

          JOBU1 is CHARACTER


          = 'Y':      U1 is computed;


          otherwise:  U1 is not computed.

JOBU2

          JOBU2 is CHARACTER


          = 'Y':      U2 is computed;


          otherwise:  U2 is not computed.

JOBV1T

          JOBV1T is CHARACTER


          = 'Y':      V1T is computed;


          otherwise:  V1T is not computed.

M

          M is INTEGER


          The number of rows in X.

P

          P is INTEGER


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

Q

          Q is INTEGER


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

X11

          X11 is REAL array, dimension (LDX11,Q)


          On entry, part of the orthogonal matrix whose CSD is desired.

LDX11

          LDX11 is INTEGER


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

X21

          X21 is REAL array, dimension (LDX21,Q)


          On entry, part of the orthogonal matrix whose CSD is desired.

LDX21

          LDX21 is INTEGER


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

THETA

          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

          U1 is REAL array, dimension (P)


          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.

LDU1

          LDU1 is INTEGER


          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=


          MAX(1,P).

U2

          U2 is REAL array, dimension (M-P)


          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal


          matrix U2.

LDU2

          LDU2 is INTEGER


          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=


          MAX(1,M-P).

V1T

          V1T is REAL array, dimension (Q)


          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal


          matrix V1**T.

LDV1T

          LDV1T is INTEGER


          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=


          MAX(1,Q).

WORK

          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

          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

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

INFO

          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.

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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. N >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQLF.

WORK

          WORK is REAL array, dimension (N)

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. N >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQRF.

WORK

          WORK is REAL array, dimension (N)

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          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

          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

          LDA is INTEGER


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

TAU

          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

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. M >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGELQF.

WORK

          WORK is REAL array, dimension (M)

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. M >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGELQF.

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= max(1,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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. N >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQLF.

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


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


          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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. N >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQRF.

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


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


          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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. M >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGERQF.

WORK

          WORK is REAL array, dimension (M)

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines the


          matrix Q. M >= K >= 0.

A

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGERQF.

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= max(1,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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          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

          N is INTEGER


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

A

          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

          LDA is INTEGER


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

TAU

          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

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK. LWORK >= max(1,N-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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left


          = 'R': apply Q or Q**T from the Right

TRANS

          TRANS is CHARACTER*1


          = 'N': apply Q  (No transpose)


          = 'T': apply Q**T (Transpose)

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQLF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                                   (N) if SIDE = 'L',


                                   (M) if SIDE = 'R'

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left


          = 'R': apply Q or Q**T from the Right

TRANS

          TRANS is CHARACTER*1


          = 'N': apply Q  (No transpose)


          = 'T': apply Q**T (Transpose)

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQRF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                                   (N) if SIDE = 'L',


                                   (M) if SIDE = 'R'

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 If 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

          VECT is CHARACTER*1


          = 'Q': apply Q or Q**T;


          = 'P': apply P or P**T.

SIDE

          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

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q  or P;


          = 'T':  Transpose, apply Q**T or P**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          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

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

ILO

          ILO is INTEGER

IHI

          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

          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

          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

          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

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left


          = 'R': apply Q or Q**T from the Right

TRANS

          TRANS is CHARACTER*1


          = 'N': apply Q  (No transpose)


          = 'T': apply Q**T (Transpose)

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGELQF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                                   (N) if SIDE = 'L',


                                   (M) if SIDE = 'R'

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGELQF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQLF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGEQRF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left


          = 'R': apply Q or Q**T from the Right

TRANS

          TRANS is CHARACTER*1


          = 'N': apply Q  (No transpose)


          = 'T': apply Q' (Transpose)

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGERQF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                                   (N) if SIDE = 'L',


                                   (M) if SIDE = 'R'

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left


          = 'R': apply Q or Q**T from the Right

TRANS

          TRANS is CHARACTER*1


          = 'N': apply Q  (No transpose)


          = 'T': apply Q**T (Transpose)

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by STZRZF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension


                                   (N) if SIDE = 'L',


                                   (M) if SIDE = 'R'

INFO

          INFO is INTEGER


          = 0: successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by SGERQF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          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

          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

          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

          LDA is INTEGER


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

TAU

          TAU is REAL array, dimension (K)


          TAU(i) must contain the scalar factor of the elementary


          reflector H(i), as returned by STZRZF.

C

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

 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

          SIDE is CHARACTER*1


          = 'L': apply Q or Q**T from the Left;


          = 'R': apply Q or Q**T from the Right.

UPLO

          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

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'T':  Transpose, apply Q**T.

M

          M is INTEGER


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

N

          N is INTEGER


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

A

          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

          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

          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

          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

          LDC is INTEGER


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

WORK

          WORK is REAL array, dimension (MAX(1,LWORK))


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

LWORK

          LWORK is INTEGER


          The dimension of the array WORK.


          If 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

          INFO is INTEGER


          = 0:  successful exit


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