CTPLQT computes a blocked LQ factorization of a complex


 'triangular-pentagonal' matrix C, which is composed of a


 triangular block A and pentagonal block B, using the compact


 WY representation for Q.

M

          M is INTEGER


          The number of rows of the matrix B, and the order of the


          triangular matrix A.


          M >= 0.

N

          N is INTEGER


          The number of columns of the matrix B.


          N >= 0.

L

          L is INTEGER


          The number of rows of the lower trapezoidal part of B.


          MIN(M,N) >= L >= 0.  See Further Details.

MB

          MB is INTEGER


          The block size to be used in the blocked QR.  M >= MB >= 1.

A

          A is COMPLEX array, dimension (LDA,M)


          On entry, the lower triangular M-by-M matrix A.


          On exit, the elements on and below the diagonal of the array


          contain the lower triangular matrix L.

LDA

          LDA is INTEGER


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

B

          B is COMPLEX array, dimension (LDB,N)


          On entry, the pentagonal M-by-N matrix B.  The first N-L columns


          are rectangular, and the last L columns are lower trapezoidal.


          On exit, B contains the pentagonal matrix V.  See Further Details.

LDB

          LDB is INTEGER


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

T

          T is COMPLEX array, dimension (LDT,N)


          The lower triangular block reflectors stored in compact form


          as a sequence of upper triangular blocks.  See Further Details.

LDT

          LDT is INTEGER


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

WORK

          WORK is COMPLEX array, dimension (MB*M)

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 input matrix C is a M-by-(M+N) matrix


               C = [ A ] [ B ]


  where A is an lower triangular M-by-M matrix, and B is M-by-N pentagonal


  matrix consisting of a M-by-(N-L) rectangular matrix B1 on left of a M-by-L


  upper trapezoidal matrix B2:


          [ B ] = [ B1 ] [ B2 ]


                   [ B1 ]  <- M-by-(N-L) rectangular


                   [ B2 ]  <-     M-by-L lower trapezoidal.


  The lower trapezoidal matrix B2 consists of the first L columns of a


  M-by-M lower triangular matrix, where 0 <= L <= MIN(M,N).  If L=0,


  B is rectangular M-by-N; if M=L=N, B is lower triangular.


  The matrix W stores the elementary reflectors H(i) in the i-th row


  above the diagonal (of A) in the M-by-(M+N) input matrix C


            [ C ] = [ A ] [ B ]


                   [ A ]  <- lower triangular M-by-M


                   [ B ]  <- M-by-N pentagonal


  so that W can be represented as


            [ W ] = [ I ] [ V ]


                   [ I ]  <- identity, M-by-M


                   [ V ]  <- M-by-N, same form as B.


  Thus, all of information needed for W is contained on exit in B, which


  we call V above.  Note that V has the same form as B; that is,


            [ V ] = [ V1 ] [ V2 ]


                   [ V1 ] <- M-by-(N-L) rectangular


                   [ V2 ] <-     M-by-L lower trapezoidal.


  The rows of V represent the vectors which define the H(i)'s.


  The number of blocks is B = ceiling(M/MB), where each


  block is of order MB except for the last block, which is of order


  IB = M - (M-1)*MB.  For each of the B blocks, a upper triangular block


  reflector factor is computed: T1, T2, ..., TB.  The MB-by-MB (and IB-by-IB


  for the last block) T's are stored in the MB-by-N matrix T as


               T = [T1 T2 ... TB].

 CTPLQT2 computes a LQ a factorization of a complex 'triangular-pentagonal'


 matrix C, which is composed of a triangular block A and pentagonal block B,


 using the compact WY representation for Q.

M

          M is INTEGER


          The total number of rows of the matrix B.


          M >= 0.

N

          N is INTEGER


          The number of columns of the matrix B, and the order of


          the triangular matrix A.


          N >= 0.

L

          L is INTEGER


          The number of rows of the lower trapezoidal part of B.


          MIN(M,N) >= L >= 0.  See Further Details.

A

          A is COMPLEX array, dimension (LDA,M)


          On entry, the lower triangular M-by-M matrix A.


          On exit, the elements on and below the diagonal of the array


          contain the lower triangular matrix L.

LDA

          LDA is INTEGER


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

B

          B is COMPLEX array, dimension (LDB,N)


          On entry, the pentagonal M-by-N matrix B.  The first N-L columns


          are rectangular, and the last L columns are lower trapezoidal.


          On exit, B contains the pentagonal matrix V.  See Further Details.

LDB

          LDB is INTEGER


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

T

          T is COMPLEX array, dimension (LDT,M)


          The N-by-N upper triangular factor T of the block reflector.


          See Further Details.

LDT

          LDT is INTEGER


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

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 input matrix C is a M-by-(M+N) matrix


               C = [ A ][ B ]


  where A is an lower triangular M-by-M matrix, and B is M-by-N pentagonal


  matrix consisting of a M-by-(N-L) rectangular matrix B1 left of a M-by-L


  upper trapezoidal matrix B2:


               B = [ B1 ][ B2 ]


                   [ B1 ]  <-     M-by-(N-L) rectangular


                   [ B2 ]  <-     M-by-L lower trapezoidal.


  The lower trapezoidal matrix B2 consists of the first L columns of a


  N-by-N lower triangular matrix, where 0 <= L <= MIN(M,N).  If L=0,


  B is rectangular M-by-N; if M=L=N, B is lower triangular.


  The matrix W stores the elementary reflectors H(i) in the i-th row


  above the diagonal (of A) in the M-by-(M+N) input matrix C


               C = [ A ][ B ]


                   [ A ]  <- lower triangular M-by-M


                   [ B ]  <- M-by-N pentagonal


  so that W can be represented as


               W = [ I ][ V ]


                   [ I ]  <- identity, M-by-M


                   [ V ]  <- M-by-N, same form as B.


  Thus, all of information needed for W is contained on exit in B, which


  we call V above.  Note that V has the same form as B; that is,


               W = [ V1 ][ V2 ]


                   [ V1 ] <-     M-by-(N-L) rectangular


                   [ V2 ] <-     M-by-L lower trapezoidal.


  The rows of V represent the vectors which define the H(i)'s.


  The (M+N)-by-(M+N) block reflector H is then given by


               H = I - W**T * T * W


  where W^H is the conjugate transpose of W and T is the upper triangular


  factor of the block reflector.

 CTPMLQT applies a complex unitary matrix Q obtained from a


 'triangular-pentagonal' complex block reflector H to a general


 complex matrix C, which consists of two blocks A and B.

SIDE

          SIDE is CHARACTER*1


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


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

TRANS

          TRANS is CHARACTER*1


          = 'N':  No transpose, apply Q;


          = 'C':  Conjugate transpose, apply Q**H.

M

          M is INTEGER


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

N

          N is INTEGER


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

K

          K is INTEGER


          The number of elementary reflectors whose product defines


          the matrix Q.

L

          L is INTEGER


          The order of the trapezoidal part of V.


          K >= L >= 0.  See Further Details.

MB

          MB is INTEGER


          The block size used for the storage of T.  K >= MB >= 1.


          This must be the same value of MB used to generate T


          in CTPLQT.

V

          V is COMPLEX array, dimension (LDV,K)


          The i-th row must contain the vector which defines the


          elementary reflector H(i), for i = 1,2,...,k, as returned by


          CTPLQT in B.  See Further Details.

LDV

          LDV is INTEGER


          The leading dimension of the array V. LDV >= K.

T

          T is COMPLEX array, dimension (LDT,K)


          The upper triangular factors of the block reflectors


          as returned by CTPLQT, stored as a MB-by-K matrix.

LDT

          LDT is INTEGER


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

A

          A is COMPLEX array, dimension


          (LDA,N) if SIDE = 'L' or


          (LDA,K) if SIDE = 'R'


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


          On exit, A is overwritten by the corresponding block of


          Q*C or Q**H*C or C*Q or C*Q**H.  See Further Details.

LDA

          LDA is INTEGER


          The leading dimension of the array A.


          If SIDE = 'L', LDA >= max(1,K);


          If SIDE = 'R', LDA >= max(1,M).

B

          B is COMPLEX array, dimension (LDB,N)


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


          On exit, B is overwritten by the corresponding block of


          Q*C or Q**H*C or C*Q or C*Q**H.  See Further Details.

LDB

          LDB is INTEGER


          The leading dimension of the array B.


          LDB >= max(1,M).

WORK

          WORK is COMPLEX array. The dimension of WORK is


           N*MB if SIDE = 'L', or  M*MB 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.

  The columns of the pentagonal matrix V contain the elementary reflectors


  H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a


  trapezoidal block V2:


        V = [V1] [V2].


  The size of the trapezoidal block V2 is determined by the parameter L,


  where 0 <= L <= K; V2 is lower trapezoidal, consisting of the first L


  rows of a K-by-K upper triangular matrix.  If L=K, V2 is lower triangular;


  if L=0, there is no trapezoidal block, hence V = V1 is rectangular.


  If SIDE = 'L':  C = [A]  where A is K-by-N,  B is M-by-N and V is K-by-M.


                      [B]


  If SIDE = 'R':  C = [A B]  where A is M-by-K, B is M-by-N and V is K-by-N.


  The complex unitary matrix Q is formed from V and T.


  If TRANS='N' and SIDE='L', C is on exit replaced with Q * C.


  If TRANS='C' and SIDE='L', C is on exit replaced with Q**H * C.


  If TRANS='N' and SIDE='R', C is on exit replaced with C * Q.


  If TRANS='C' and SIDE='R', C is on exit replaced with C * Q**H.

 DBBCSD 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 DORCSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (Q)


          When DBBCSD converges, B11D contains the cosines of THETA(1),


          ..., THETA(Q). If DBBCSD fails to converge, then B11D


          contains the diagonal of the partially reduced top-left


          block.

B11E

          B11E is DOUBLE PRECISION array, dimension (Q-1)


          When DBBCSD converges, B11E contains zeros. If DBBCSD fails


          to converge, then B11E contains the superdiagonal of the


          partially reduced top-left block.

B12D

          B12D is DOUBLE PRECISION array, dimension (Q)


          When DBBCSD converges, B12D contains the negative sines of


          THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then


          B12D contains the diagonal of the partially reduced top-right


          block.

B12E

          B12E is DOUBLE PRECISION array, dimension (Q-1)


          When DBBCSD converges, B12E contains zeros. If DBBCSD fails


          to converge, then B12E contains the subdiagonal of the


          partially reduced top-right block.

B21D

          B21D is DOUBLE PRECISION  array, dimension (Q)


          When DBBCSD converges, B21D contains the negative sines of


          THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then


          B21D contains the diagonal of the partially reduced bottom-left


          block.

B21E

          B21E is DOUBLE PRECISION  array, dimension (Q-1)


          When DBBCSD converges, B21E contains zeros. If DBBCSD fails


          to converge, then B21E contains the subdiagonal of the


          partially reduced bottom-left block.

B22D

          B22D is DOUBLE PRECISION  array, dimension (Q)


          When DBBCSD converges, B22D contains the negative sines of


          THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then


          B22D contains the diagonal of the partially reduced bottom-right


          block.

B22E

          B22E is DOUBLE PRECISION  array, dimension (Q-1)


          When DBBCSD converges, B22E contains zeros. If DBBCSD fails


          to converge, then B22E contains the subdiagonal of the


          partially reduced bottom-right block.

WORK

          WORK is DOUBLE PRECISION 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 DBBCSD did not converge, INFO specifies the number


                of nonzero entries in PHI, and B11D, B11E, etc.,


                contain the partially reduced matrix.

  TOLMUL  DOUBLE PRECISION, 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.

 DGETSQRHRT computes a NB2-sized column blocked QR-factorization


 of a real M-by-N matrix A with M >= N,


    A = Q * R.


 The routine uses internally a NB1-sized column blocked and MB1-sized


 row blocked TSQR-factorization and perfors the reconstruction


 of the Householder vectors from the TSQR output. The routine also


 converts the R_tsqr factor from the TSQR-factorization output into


 the R factor that corresponds to the Householder QR-factorization,


    A = Q_tsqr * R_tsqr = Q * R.


 The output Q and R factors are stored in the same format as in DGEQRT


 (Q is in blocked compact WY-representation). See the documentation


 of DGEQRT for more details on the format.

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. M >= N >= 0.

MB1

          MB1 is INTEGER


          The row block size to be used in the blocked TSQR.


          MB1 > N.

NB1

          NB1 is INTEGER


          The column block size to be used in the blocked TSQR.


          N >= NB1 >= 1.

NB2

          NB2 is INTEGER


          The block size to be used in the blocked QR that is


          output. NB2 >= 1.

A

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


          On entry: an M-by-N matrix A.


          On exit:


           a) the elements on and above the diagonal


              of the array contain the N-by-N upper-triangular


              matrix R corresponding to the Householder QR;


           b) the elements below the diagonal represent Q by


              the columns of blocked V (compact WY-representation).

LDA

          LDA is INTEGER


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

T

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


          The upper triangular block reflectors stored in compact form


          as a sequence of upper triangular blocks.

LDT

          LDT is INTEGER


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

WORK

          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))


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

LWORK

          The dimension of the array WORK.


          LWORK >= MAX( LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) ),


          where


             NUM_ALL_ROW_BLOCKS = CEIL((M-N)/(MB1-N)),


             NB1LOCAL = MIN(NB1,N).


             LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL,


             LW1 = NB1LOCAL * N,


             LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) ),


          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.

 November 2020, Igor Kozachenko,


                Computer Science Division,


                University of California, Berkeley

 DGGHD3 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 DGGHD3 reduces the original


 problem to generalized Hessenberg form.


 This is a blocked variant of DGGHRD, 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 DGGBAL; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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).

 DGGHRD 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 DGGHRD 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 DGGBAL; 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.)

 DGGQRF 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(N,M))


          The scalar factors of the elementary reflectors which


          represent the orthogonal matrix Q (see Further Details).

B

          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(N,P))


          The scalar factors of the elementary reflectors which


          represent the orthogonal matrix Z (see Further Details).

WORK

          WORK is DOUBLE PRECISION 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 DORMQR.


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


  To use Q to update another matrix, use LAPACK subroutine DORMQR.


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


  To use Z to update another matrix, use LAPACK subroutine DORMRQ.

 DGGRQF 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(M,N))


          The scalar factors of the elementary reflectors which


          represent the orthogonal matrix Q (see Further Details).

B

          B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(P,N))


          The scalar factors of the elementary reflectors which


          represent the orthogonal matrix Z (see Further Details).

WORK

          WORK is DOUBLE PRECISION 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 DORMRQ.


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


  To use Q to update another matrix, use LAPACK subroutine DORMRQ.


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


  To use Z to update another matrix, use LAPACK subroutine DORMQR.

 DGGSVP3 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


 DGGSVD3.

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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION

TOLB

          TOLB is DOUBLE PRECISION


          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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)

WORK

          WORK is DOUBLE PRECISION 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 DGEQP3 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.


  DGGSVP3 replaces the deprecated subroutine DGGSVP.

 DGSVJ0 is called from DGESVJ as a pre-processor and that is its main


 purpose. It applies Jacobi rotations in the same way as DGESVJ 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION


          EPS = DLAMCH('Epsilon')

SFMIN

          SFMIN is DOUBLE PRECISION


          SFMIN = DLAMCH('Safe Minimum')

TOL

          TOL is DOUBLE PRECISION


          TOL is the threshold for Jacobi rotations. For a pair


          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is


          applied only if DABS(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 DOUBLE PRECISION 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.

DGSVJ0 is used just to enable DGESVJ 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.

 DGSVJ1 is called from DGESVJ as a pre-processor and that is its main


 purpose. It applies Jacobi rotations in the same way as DGESVJ 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


 ~~~~~~~~~~~~~~~


 DGSVJ1 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION


          EPS = DLAMCH('Epsilon')

SFMIN

          SFMIN is DOUBLE PRECISION


          SFMIN = DLAMCH('Safe Minimum')

TOL

          TOL is DOUBLE PRECISION


          TOL is the threshold for Jacobi rotations. For a pair


          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is


          applied only if DABS(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 DOUBLE PRECISION 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)

 DHSEIN 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 DHSEQR; 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 DHSEIN 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, DHSEIN 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)

WI

          WI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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|.

    DHSEQR 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 DGEBAL, and then passed to ZGEHRD


           when the matrix output by DGEBAL is reduced to Hessenberg


           form. Otherwise ILO and IHI should be set to 1 and N


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


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

H

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


           On entry, the upper Hessenberg matrix H.


           On exit, if INFO = 0 and 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 DHSEQR, 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 DOUBLE PRECISION array, dimension (N)

WI

          WI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DORGHR


           after the call to DGEHRD 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 DOUBLE PRECISION 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 DHSEQR does a workspace query.


           In this case, DHSEQR 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, DHSEQR 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,'DHSEQR',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 DLAHQR vs DLAQR0 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


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

    DLA_LIN_BERR computes component-wise 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 component-wise 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 DOUBLE PRECISION array, dimension (N,NRHS)


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


     error formula above.

AYB

          AYB is DOUBLE PRECISION 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 dla_gerfsx_extended.f).

BERR

          BERR is DOUBLE PRECISION array, dimension (NRHS)


     The component-wise relative backward error from the formula above.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

    DLA_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 DOUBLE PRECISION array, dimension (N)


            The first part of the doubled-single accumulation vector.

Y

          Y is DOUBLE PRECISION array, dimension (N)


            The second part of the doubled-single accumulation vector.

W

          W is DOUBLE PRECISION array, dimension (N)


            The vector to be added.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 DLALS0 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION


         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 DOUBLE PRECISION


         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 DOUBLE PRECISION 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

 DLALSA 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, DLALSA applies the inverse of the left singular vector


 matrix of an upper bidiagonal matrix to the right hand side; and if


 ICOMPQ = 1, DLALSA applies the right singular vector matrix to the


 right hand side. The singular vector matrices were generated in


 compact form by DLALSA.

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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension ( LDU, NLVL ).


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

DIFR

          DIFR is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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

 DLALSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N-1)


         Contains the super-diagonal entries of the bidiagonal matrix.


         On exit, E has been destroyed.

B

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


         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 DOUBLE PRECISION 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

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

DLANSF

    DLANSF = ( 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 DLANSF 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, DLANSF is


          set to zero.

A

          A is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (MAX(1,LWORK)),


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


          WORK is not referenced.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

  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

 DLARSCL2 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_dge_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 DOUBLE PRECISION array, dimension (M)


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

X

          X is DOUBLE PRECISION 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.

 DLARZ 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 DTZRZF.

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 DOUBLE PRECISION array, dimension (1+(L-1)*abs(INCV))


          The vector v in the representation of H as returned by


          DTZRZF. V is not used if TAU = 0.

INCV

          INCV is INTEGER


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

TAU

          TAU is DOUBLE PRECISION


          The value tau in the representation of H.

C

          C is DOUBLE PRECISION array, dimension (LDC,N)


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


          On exit, C is overwritten by the matrix H * C if SIDE = 'L',


          or C * H if SIDE = 'R'.

LDC

          LDC is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension


                         (N) if SIDE = 'L'


                      or (M) if SIDE = 'R'

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

 DLARZB 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDT,K)


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


          block reflector.

LDT

          LDT is INTEGER


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

C

          C is DOUBLE PRECISION array, dimension (LDC,N)


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


          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.

LDC

          LDC is INTEGER


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

WORK

          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)

LDWORK

          LDWORK is INTEGER


          The leading dimension of the array WORK.


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


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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

 DLARZT 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 DOUBLE PRECISION array, dimension


                               (LDV,K) if STOREV = 'C'


                               (LDV,N) if STOREV = 'R'


          The matrix V. See further details.

LDV

          LDV is INTEGER


          The leading dimension of the array V.


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

TAU

          TAU is DOUBLE PRECISION array, dimension (K)


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


          reflector H(i).

T

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


          The k by k triangular factor T of the block reflector.


          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is


          lower triangular. The rest of the array is not used.

LDT

          LDT is INTEGER


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

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

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 )

 DLASCL2 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_dge_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 DOUBLE PRECISION array, length M


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

X

          X is DOUBLE PRECISION 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.

 DLATRZ 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (M)


          The scalar factors of the elementary reflectors.

WORK

          WORK is DOUBLE PRECISION 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 ).

 DOPGTR 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


 DSPTRD 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 DSPTRD;


          = 'L': Lower triangular packed storage used in previous


                 call to DSPTRD.

N

          N is INTEGER


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

AP

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


          The vectors which define the elementary reflectors, as


          returned by DSPTRD.

TAU

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


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


          reflector H(i), as returned by DSPTRD.

Q

          Q is DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DOPMTR 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 DSPTRD 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 DSPTRD;


          = 'L': Lower triangular packed storage used in previous


                 call to DSPTRD.

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 DOUBLE PRECISION 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 DSPTRD.  AP is modified by the routine but


          restored on exit.

TAU

          TAU is DOUBLE PRECISION 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 DSPTRD.

C

          C is DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DORBDB 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 DORCSD


 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (P)


          The scalar factors of the elementary reflectors that define


          P1.

TAUP2

          TAUP2 is DOUBLE PRECISION array, dimension (M-P)


          The scalar factors of the elementary reflectors that define


          P2.

TAUQ1

          TAUQ1 is DOUBLE PRECISION array, dimension (Q)


          The scalar factors of the elementary reflectors that define


          Q1.

TAUQ2

          TAUQ2 is DOUBLE PRECISION array, dimension (M-Q)


          The scalar factors of the elementary reflectors that define


          Q2.

WORK

          WORK is DOUBLE PRECISION 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 DORCSD for details.


  P1, P2, Q1, and Q2 are represented as products of elementary


  reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2


  using DORGQR and DORGLQ.

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

 DORBDB1 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 DORBDB2, DORBDB3, and DORBDB4 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (Q)


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


           THETA and PHI. See Further Details.

PHI

          PHI is DOUBLE PRECISION array, dimension (Q-1)


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


           THETA and PHI. See Further Details.

TAUP1

          TAUP1 is DOUBLE PRECISION array, dimension (P)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is DOUBLE PRECISION array, dimension (M-P)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is DOUBLE PRECISION array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

WORK

          WORK is DOUBLE PRECISION 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 DORCSD for details.


  P1, P2, and Q1 are represented as products of elementary reflectors.


  See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR


  and DORGLQ.

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

 DORBDB2 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 DORBDB1, DORBDB3, and DORBDB4 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (Q)


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


           THETA and PHI. See Further Details.

PHI

          PHI is DOUBLE PRECISION array, dimension (Q-1)


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


           THETA and PHI. See Further Details.

TAUP1

          TAUP1 is DOUBLE PRECISION array, dimension (P-1)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is DOUBLE PRECISION array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is DOUBLE PRECISION array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

WORK

          WORK is DOUBLE PRECISION 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 DORCSD for details.


  P1, P2, and Q1 are represented as products of elementary reflectors.


  See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR


  and DORGLQ.

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

 DORBDB3 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 DORBDB1, DORBDB2, and DORBDB4 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (Q)


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


           THETA and PHI. See Further Details.

PHI

          PHI is DOUBLE PRECISION array, dimension (Q-1)


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


           THETA and PHI. See Further Details.

TAUP1

          TAUP1 is DOUBLE PRECISION array, dimension (P)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is DOUBLE PRECISION array, dimension (M-P)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is DOUBLE PRECISION array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

WORK

          WORK is DOUBLE PRECISION 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 DORCSD for details.


  P1, P2, and Q1 are represented as products of elementary reflectors.


  See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR


  and DORGLQ.

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

 DORBDB4 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 DORBDB1, DORBDB2, and DORBDB3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (Q)


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


           THETA and PHI. See Further Details.

PHI

          PHI is DOUBLE PRECISION array, dimension (Q-1)


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


           THETA and PHI. See Further Details.

TAUP1

          TAUP1 is DOUBLE PRECISION array, dimension (M-Q)


           The scalar factors of the elementary reflectors that define


           P1.

TAUP2

          TAUP2 is DOUBLE PRECISION array, dimension (M-Q)


           The scalar factors of the elementary reflectors that define


           P2.

TAUQ1

          TAUQ1 is DOUBLE PRECISION array, dimension (Q)


           The scalar factors of the elementary reflectors that define


           Q1.

PHANTOM

          PHANTOM is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DORCSD for details.


  P1, P2, and Q1 are represented as products of elementary reflectors.


  See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR


  and DORGLQ.

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

 DORBDB5 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DORBDB6 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DORCSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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:  DBBCSD 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.

 DORCSD2BY1 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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:  DBBCSD 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.

 DORG2L 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 DGEQLF.

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 DOUBLE PRECISION 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 DGEQLF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGEQLF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORG2R 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 DGEQRF.

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 DOUBLE PRECISION 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 DGEQRF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGEQRF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGHR generates a real orthogonal matrix Q which is defined as the


 product of IHI-ILO elementary reflectors of order N, as returned by


 DGEHRD:


 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 DGEHRD. 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 DOUBLE PRECISION array, dimension (LDA,N)


          On entry, the vectors which define the elementary reflectors,


          as returned by DGEHRD.


          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 DOUBLE PRECISION array, dimension (N-1)


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


          reflector H(i), as returned by DGEHRD.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGL2 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 DGELQF.

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 DOUBLE PRECISION 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 DGELQF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGELQF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGLQ 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 DGELQF.

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 DOUBLE PRECISION 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 DGELQF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGELQF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGQL 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 DGEQLF.

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 DOUBLE PRECISION 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 DGEQLF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGEQLF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGQR 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 DGEQRF.

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 DOUBLE PRECISION 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 DGEQRF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGEQRF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGR2 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 DGERQF.

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 DOUBLE PRECISION 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 DGERQF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGERQF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGRQ 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 DGERQF.

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 DOUBLE PRECISION 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 DGERQF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGERQF.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGTR generates a real orthogonal matrix Q which is defined as the


 product of n-1 elementary reflectors of order N, as returned by


 DSYTRD:


 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 DSYTRD;


          = 'L': Lower triangle of A contains elementary reflectors


                 from DSYTRD.

N

          N is INTEGER


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

A

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


          On entry, the vectors which define the elementary reflectors,


          as returned by DSYTRD.


          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 DOUBLE PRECISION array, dimension (N-1)


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


          reflector H(i), as returned by DSYTRD.

WORK

          WORK is DOUBLE PRECISION 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.

 DORGTSQR generates an M-by-N real matrix Q_out with orthonormal columns,


 which are the first N columns of a product of real orthogonal


 matrices of order M which are returned by DLATSQR


      Q_out = first_N_columns_of( Q(1)_in * Q(2)_in * ... * Q(k)_in ).


 See the documentation for DLATSQR.

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. M >= N >= 0.

MB

          MB is INTEGER


          The row block size used by DLATSQR to return


          arrays A and T. MB > N.


          (Note that if MB > M, then M is used instead of MB


          as the row block size).

NB

          NB is INTEGER


          The column block size used by DLATSQR to return


          arrays A and T. NB >= 1.


          (Note that if NB > N, then N is used instead of NB


          as the column block size).

A

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


          On entry:


             The elements on and above the diagonal are not accessed.


             The elements below the diagonal represent the unit


             lower-trapezoidal blocked matrix V computed by DLATSQR


             that defines the input matrices Q_in(k) (ones on the


             diagonal are not stored) (same format as the output A


             below the diagonal in DLATSQR).


          On exit:


             The array A contains an M-by-N orthonormal matrix Q_out,


             i.e the columns of A are orthogonal unit vectors.

LDA

          LDA is INTEGER


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

T

          T is DOUBLE PRECISION array,


          dimension (LDT, N * NIRB)


          where NIRB = Number_of_input_row_blocks


                     = MAX( 1, CEIL((M-N)/(MB-N)) )


          Let NICB = Number_of_input_col_blocks


                   = CEIL(N/NB)


          The upper-triangular block reflectors used to define the


          input matrices Q_in(k), k=(1:NIRB*NICB). The block


          reflectors are stored in compact form in NIRB block


          reflector sequences. Each of NIRB block reflector sequences


          is stored in a larger NB-by-N column block of T and consists


          of NICB smaller NB-by-NB upper-triangular column blocks.


          (same format as the output T in DLATSQR).

LDT

          LDT is INTEGER


          The leading dimension of the array T.


          LDT >= max(1,min(NB1,N)).

WORK

          (workspace) DOUBLE PRECISION array, dimension (MAX(2,LWORK))


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

LWORK

          The dimension of the array WORK.  LWORK >= (M+NB)*N.


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


          The routine only calculates the optimal size of the WORK


          array, returns this value as the first entry of the WORK


          array, and no error message related to LWORK is issued


          by XERBLA.

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.

 November 2019, Igor Kozachenko,


                Computer Science Division,


                University of California, Berkeley

 DORGTSQR_ROW generates an M-by-N real matrix Q_out with


 orthonormal columns from the output of DLATSQR. These N orthonormal


 columns are the first N columns of a product of complex unitary


 matrices Q(k)_in of order M, which are returned by DLATSQR in


 a special format.


      Q_out = first_N_columns_of( Q(1)_in * Q(2)_in * ... * Q(k)_in ).


 The input matrices Q(k)_in are stored in row and column blocks in A.


 See the documentation of DLATSQR for more details on the format of


 Q(k)_in, where each Q(k)_in is represented by block Householder


 transformations. This routine calls an auxiliary routine DLARFB_GETT,


 where the computation is performed on each individual block. The


 algorithm first sweeps NB-sized column blocks from the right to left


 starting in the bottom row block and continues to the top row block


 (hence _ROW in the routine name). This sweep is in reverse order of


 the order in which DLATSQR generates the output blocks.

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. M >= N >= 0.

MB

          MB is INTEGER


          The row block size used by DLATSQR to return


          arrays A and T. MB > N.


          (Note that if MB > M, then M is used instead of MB


          as the row block size).

NB

          NB is INTEGER


          The column block size used by DLATSQR to return


          arrays A and T. NB >= 1.


          (Note that if NB > N, then N is used instead of NB


          as the column block size).

A

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


          On entry:


             The elements on and above the diagonal are not used as


             input. The elements below the diagonal represent the unit


             lower-trapezoidal blocked matrix V computed by DLATSQR


             that defines the input matrices Q_in(k) (ones on the


             diagonal are not stored). See DLATSQR for more details.


          On exit:


             The array A contains an M-by-N orthonormal matrix Q_out,


             i.e the columns of A are orthogonal unit vectors.

LDA

          LDA is INTEGER


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

T

          T is DOUBLE PRECISION array,


          dimension (LDT, N * NIRB)


          where NIRB = Number_of_input_row_blocks


                     = MAX( 1, CEIL((M-N)/(MB-N)) )


          Let NICB = Number_of_input_col_blocks


                   = CEIL(N/NB)


          The upper-triangular block reflectors used to define the


          input matrices Q_in(k), k=(1:NIRB*NICB). The block


          reflectors are stored in compact form in NIRB block


          reflector sequences. Each of the NIRB block reflector


          sequences is stored in a larger NB-by-N column block of T


          and consists of NICB smaller NB-by-NB upper-triangular


          column blocks. See DLATSQR for more details on the format


          of T.

LDT

          LDT is INTEGER


          The leading dimension of the array T.


          LDT >= max(1,min(NB,N)).

WORK

          (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))


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

LWORK

          The dimension of the array WORK.


          LWORK >= NBLOCAL * MAX(NBLOCAL,(N-NBLOCAL)),


          where NBLOCAL=MIN(NB,N).


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


          The routine only calculates the optimal size of the WORK


          array, returns this value as the first entry of the WORK


          array, and no error message related to LWORK is issued


          by XERBLA.

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.

 November 2020, Igor Kozachenko,


                Computer Science Division,


                University of California, Berkeley

  DORHR_COL takes an M-by-N real matrix Q_in with orthonormal columns


  as input, stored in A, and performs Householder Reconstruction (HR),


  i.e. reconstructs Householder vectors V(i) implicitly representing


  another M-by-N matrix Q_out, with the property that Q_in = Q_out*S,


  where S is an N-by-N diagonal matrix with diagonal entries


  equal to +1 or -1. The Householder vectors (columns V(i) of V) are


  stored in A on output, and the diagonal entries of S are stored in D.


  Block reflectors are also returned in T


  (same output format as DGEQRT).

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. M >= N >= 0.

NB

          NB is INTEGER


          The column block size to be used in the reconstruction


          of Householder column vector blocks in the array A and


          corresponding block reflectors in the array T. NB >= 1.


          (Note that if NB > N, then N is used instead of NB


          as the column block size.)

A

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


          On entry:


             The array A contains an M-by-N orthonormal matrix Q_in,


             i.e the columns of A are orthogonal unit vectors.


          On exit:


             The elements below the diagonal of A represent the unit


             lower-trapezoidal matrix V of Householder column vectors


             V(i). The unit diagonal entries of V are not stored


             (same format as the output below the diagonal in A from


             DGEQRT). The matrix T and the matrix V stored on output


             in A implicitly define Q_out.


             The elements above the diagonal contain the factor U


             of the 'modified' LU-decomposition:


                Q_in - ( S ) = V * U


                       ( 0 )


             where 0 is a (M-N)-by-(M-N) zero matrix.

LDA

          LDA is INTEGER


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

T

          T is DOUBLE PRECISION array,


          dimension (LDT, N)


          Let NOCB = Number_of_output_col_blocks


                   = CEIL(N/NB)


          On exit, T(1:NB, 1:N) contains NOCB upper-triangular


          block reflectors used to define Q_out stored in compact


          form as a sequence of upper-triangular NB-by-NB column


          blocks (same format as the output T in DGEQRT).


          The matrix T and the matrix V stored on output in A


          implicitly define Q_out. NOTE: The lower triangles


          below the upper-triangular blocks will be filled with


          zeros. See Further Details.

LDT

          LDT is INTEGER


          The leading dimension of the array T.


          LDT >= max(1,min(NB,N)).

D

          D is DOUBLE PRECISION array, dimension min(M,N).


          The elements can be only plus or minus one.


          D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where


          1 <= i <= min(M,N), and Q_in_i is Q_in after performing


          i-1 steps of “modified” Gaussian elimination.


          See Further Details.

INFO

          INFO is INTEGER


          = 0:  successful exit


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

 The computed M-by-M orthogonal factor Q_out is defined implicitly as


 a product of orthogonal matrices Q_out(i). Each Q_out(i) is stored in


 the compact WY-representation format in the corresponding blocks of


 matrices V (stored in A) and T.


 The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N


 matrix A contains the column vectors V(i) in NB-size column


 blocks VB(j). For example, VB(1) contains the columns


 V(1), V(2), ... V(NB). NOTE: The unit entries on


 the diagonal of Y are not stored in A.


 The number of column blocks is


     NOCB = Number_of_output_col_blocks = CEIL(N/NB)


 where each block is of order NB except for the last block, which


 is of order LAST_NB = N - (NOCB-1)*NB.


 For example, if M=6,  N=5 and NB=2, the matrix V is


     V = (    VB(1),   VB(2), VB(3) ) =


       = (   1                      )


         ( v21    1                 )


         ( v31  v32    1            )


         ( v41  v42  v43   1        )


         ( v51  v52  v53  v54    1  )


         ( v61  v62  v63  v54   v65 )


 For each of the column blocks VB(i), an upper-triangular block


 reflector TB(i) is computed. These blocks are stored as


 a sequence of upper-triangular column blocks in the NB-by-N


 matrix T. The size of each TB(i) block is NB-by-NB, except


 for the last block, whose size is LAST_NB-by-LAST_NB.


 For example, if M=6,  N=5 and NB=2, the matrix T is


     T  = (    TB(1),    TB(2), TB(3) ) =


        = ( t11  t12  t13  t14   t15  )


          (      t22       t24        )


 The M-by-M factor Q_out is given as a product of NOCB


 orthogonal M-by-M matrices Q_out(i).


     Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB),


 where each matrix Q_out(i) is given by the WY-representation


 using corresponding blocks from the matrices V and T:


     Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T,


 where I is the identity matrix. Here is the formula with matrix


 dimensions:


  Q(i){M-by-M} = I{M-by-M} -


    VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M},


 where INB = NB, except for the last block NOCB


 for which INB=LAST_NB.


 =====


 NOTE:


 =====


 If Q_in is the result of doing a QR factorization


 B = Q_in * R_in, then:


 B = (Q_out*S) * R_in = Q_out * (S * R_in) = Q_out * R_out.


 So if one wants to interpret Q_out as the result


 of the QR factorization of B, then the corresponding R_out


 should be equal to R_out = S * R_in, i.e. some rows of R_in


 should be multiplied by -1.


 For the details of the algorithm, see [1].


 [1] 'Reconstructing Householder vectors from tall-skinny QR',


     G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,


     E. Solomonik, J. Parallel Distrib. Comput.,


     vol. 85, pp. 3-31, 2015.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

 November   2019, Igor Kozachenko,


            Computer Science Division,


            University of California, Berkeley

 DORM2L 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 DGEQLF. 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 DOUBLE PRECISION 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


          DGEQLF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGEQLF.

C

          C is DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DORM2R 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 DGEQRF. 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 DOUBLE PRECISION 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


          DGEQRF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGEQRF.

C

          C is DOUBLE PRECISION 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 DOUBLE PRECISION 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', DORMBR 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', DORMBR 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 DGEBRD 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 DGEBRD.


          If VECT = 'P', the number of rows in the original


          matrix reduced by DGEBRD.


          K >= 0.

A

          A is DOUBLE PRECISION 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 DGEBRD.

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 DOUBLE PRECISION 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 DGEBRD in the array argument TAUQ or TAUP.

C

          C is DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DORMHR 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 DGEHRD:


 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 DGEHRD. 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 DOUBLE PRECISION array, dimension


                               (LDA,M) if SIDE = 'L'


                               (LDA,N) if SIDE = 'R'


          The vectors which define the elementary reflectors, as


          returned by DGEHRD.

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 DOUBLE PRECISION 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 DGEHRD.

C

          C is DOUBLE PRECISION 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 DOUBLE PRECISION 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.

 DORML2 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 DGELQF. 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 DOUBLE PRECISION 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


          DGELQF 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 DOUBLE PRECISION array, dimension (K)


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


          reflector H(i), as returned by DGELQF.

C

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