1DGGRQF ‐ a generalized RQ factorization of an M‐by‐N matrix A and
2a P‐by‐N matrix B SUBROUTINE DGGRQF( M, P, N, A,  LDA,  TAUA,  B,
3LDB, TAUB, WORK, LWORK, INFO )
4    INTEGER INFO, LDA, LDB, LWORK, M, N, P
5    DOUBLE PRECISION A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( *
6), WORK( * ) DGGRQF computes a generalized RQ factorization of an
7M‐by‐N matrix A and a P‐by‐N matrix B:
8
9            A = R*Q,        B = Z*T*Q,
10
11where  Q is an N‐by‐N orthogonal matrix, Z is a P‐by‐P orthogonal
12matrix, and R and T assume one of the forms:
13
14if M <= N,  R = ( 0  R12 ) M,   or if M > N,  R = ( R11 ) M‐N,
15                 N‐M  M                           ( R21 ) N
16                                                     N
17
18where R12 or R21 is upper triangular, and
19
20if P >= N,  T = ( T11 ) N  ,   or if P < N,  T = ( T11  T12 ) P,
21                (  0  ) P‐N                         P   N‐P
22                   N
23
24where T11 is upper triangular.
25
26In particular, if B is square and nonsingular, the GRQ factoriza‐
27tion  of  A  and  B  implicitly  gives  the  RQ  factorization of
28A*inv(B):
29
30             A*inv(B) = (R*inv(T))*Z'
31
32where inv(B) denotes the inverse of the matrix B, and Z'  denotes
33the transpose of the matrix Z.
34
35M       (input) INTEGER The number of rows of the matrix A.  M >=
360.  P       (input) INTEGER The number of rows of the  matrix  B.
37P >= 0.  N       (input) INTEGER The number of columns of the ma‐
38trices A and B. N >= 0.  A       (input/output) DOUBLE  PRECISION
39array, dimension (LDA,N) On entry, the M‐by‐N matrix A.  On exit,
40if M <= N, the upper triangle of the subarray A(1:M,N‐M+1:N) con‐
41tains  the  M‐by‐M  upper triangular matrix R; if M > N, the ele‐
42ments on and above the (M‐N)‐th subdiagonal  contain  the  M‐by‐N
43upper  trapezoidal matrix R; the remaining elements, with the ar‐
44ray TAUA, represent the orthogonal matrix Q as a product of  ele‐
45mentary  reflectors (see Further Details).  LDA     (input) INTE‐
46GER The leading dimension of the array A. LDA >= max(1,M).   TAUA
47(output)  DOUBLE PRECISION array, dimension (min(M,N)) The scalar
48factors of the elementary reflectors which represent the orthogo‐
49nal  matrix Q (see Further Details).  B       (input/output) DOU‐
50BLE PRECISION array, dimension (LDB,N) On entry, the  P‐by‐N  ma‐
51trix  B.   On exit, the elements on and above the diagonal of the
52array contain the min(P,N)‐by‐N upper trapezoidal matrix T (T  is
53upper  triangular  if  P  >= N); the elements below the diagonal,
54with the array TAUB, represent the orthogonal matrix Z as a prod‐
55uct of elementary reflectors (see Further Details).  LDB     (in‐
56put) INTEGER The  leading  dimension  of  the  array  B.  LDB  >=
57max(1,P).   TAUB     (output)  DOUBLE  PRECISION array, dimension
58(min(P,N)) The scalar factors of the elementary reflectors  which
59represent  the  orthogonal  matrix Z (see Further Details).  WORK
60(workspace/output)    DOUBLE    PRECISION    array,     dimension
61(MAX(1,LWORK))  On exit, if INFO = 0, WORK(1) returns the optimal
62LWORK.  LWORK   (input) INTEGER The dimension of the array  WORK.
63LWORK   >=   max(1,N,M,P).   For  optimum  performance  LWORK  >=
64max(N,M,P)*max(NB1,NB2,NB3), where NB1 is the  optimal  blocksize
65for  the RQ factorization of an M‐by‐N matrix, NB2 is the optimal
66blocksize for the QR factorization of a P‐by‐N matrix, and NB3 is
67the optimal blocksize for a call of DORMRQ.
68
69If LWORK = ‐1, then a workspace query is assumed; the routine on‐
70ly calculates the optimal size of the WORK  array,  returns  this
71value  as the first entry of the WORK array, and no error message
72related to LWORK is issued by XERBLA.  INFO    (output) INTEGER =
730:  successful exit
74<  0:   if INF0= ‐i, the i‐th argument had an illegal value.  The
75matrix Q is represented as a product of elementary reflectors
76
77   Q = H(1) H(2) . . . H(k), where k = min(m,n).
78
79Each H(i) has the form
80
81   H(i) = I ‐ taua * v * v'
82
83where taua is a real scalar, and v is a real vector with
84v(n‐k+i+1:n) = 0 and v(n‐k+i) = 1; v(1:n‐k+i‐1) is stored on exit
85in A(m‐k+i,1:n‐k+i‐1), and taua in TAUA(i).
86To form Q explicitly, use LAPACK subroutine DORGRQ.
87To use Q to update another matrix, use LAPACK subroutine DORMRQ.
88
89The matrix Z is represented as a product of elementary reflectors
90
91   Z = H(1) H(2) . . . H(k), where k = min(p,n).
92
93Each H(i) has the form
94
95   H(i) = I ‐ taub * v * v'
96
97where taub is a real scalar, and v is a real vector with
98v(1:i‐1)  =  0  and  v(i)  =  1;  v(i+1:p)  is  stored on exit in
99B(i+1:p,i), and taub in TAUB(i).
100To form Z explicitly, use LAPACK subroutine DORGQR.
101To use Z to update another matrix, use LAPACK subroutine DORMQR.
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132