1ZGGRQF ‐ a generalized RQ factorization of an M‐by‐N matrix A and
2a P‐by‐N matrix B SUBROUTINE ZGGRQF( M, P, N, A,  LDA,  TAUA,  B,
3LDB, TAUB, WORK, LWORK, INFO )
4    INTEGER INFO, LDA, LDB, LWORK, M, N, P
5    COMPLEX*16  A(  LDA,  * ), B( LDB, * ), TAUA( * ), TAUB( * ),
6WORK( * ) ZGGRQF computes a generalized RQ factorization of an M‐
7by‐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 unitary matrix, Z is a P‐by‐P unitary ma‐
12trix, 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 conjugate 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) COMPLEX*16 array,
39dimension  (LDA,N)  On entry, the M‐by‐N matrix A.  On exit, if M
40<= N, the upper triangle of the subarray A(1:M,N‐M+1:N)  contains
41the  M‐by‐M  upper triangular matrix R; if M > N, the elements on
42and above the  (M‐N)‐th  subdiagonal  contain  the  M‐by‐N  upper
43trapezoidal  matrix  R;  the  remaining  elements, with the array
44TAUA, represent the unitary matrix Q as a product  of  elementary
45reflectors  (see  Further  Details).  LDA     (input) INTEGER The
46leading dimension of the array A. LDA >= max(1,M).  TAUA    (out‐
47put) COMPLEX*16 array, dimension (min(M,N)) The scalar factors of
48the elementary reflectors which represent the  unitary  matrix  Q
49(see  Further Details).  B       (input/output) COMPLEX*16 array,
50dimension (LDB,N) On entry, the P‐by‐N matrix B.   On  exit,  the
51elements  on  and  above  the  diagonal  of the array contain the
52min(P,N)‐by‐N upper trapezoidal matrix T (T is  upper  triangular
53if P >= N); the elements below the diagonal, with the array TAUB,
54represent the unitary matrix Z as a product of elementary reflec‐
55tors  (see Further Details).  LDB     (input) INTEGER The leading
56dimension of the array B. LDB >= max(1,P).  TAUB    (output) COM‐
57PLEX*16 array, dimension (min(P,N)) The scalar factors of the el‐
58ementary reflectors which represent the  unitary  matrix  Z  (see
59Further  Details).   WORK    (workspace/output) COMPLEX*16 array,
60dimension (MAX(1,LWORK)) On exit, if INFO =  0,  WORK(1)  returns
61the  optimal LWORK.  LWORK   (input) INTEGER The dimension of the
62array WORK. LWORK >= max(1,N,M,P).  For optimum performance LWORK
63>=  max(N,M,P)*max(NB1,NB2,NB3),  where NB1 is the optimal block‐
64size for the RQ factorization of an M‐by‐N matrix, NB2 is the op‐
65timal  blocksize for the QR factorization of a P‐by‐N matrix, and
66NB3 is the optimal blocksize for a call of ZUNMRQ.
67
68If LWORK = ‐1, then a workspace query is assumed; the routine on‐
69ly  calculates  the  optimal size of the WORK array, returns this
70value as the first entry of the WORK array, and no error  message
71related to LWORK is issued by XERBLA.  INFO    (output) INTEGER =
720:  successful exit
73< 0:  if INFO=‐i, the i‐th argument had an  illegal  value.   The
74matrix Q is represented as a product of elementary reflectors
75
76   Q = H(1) H(2) . . . H(k), where k = min(m,n).
77
78Each H(i) has the form
79
80   H(i) = I ‐ taua * v * v'
81
82where  taua  is  a complex scalar, and v is a complex vector with
83v(n‐k+i+1:n) = 0 and v(n‐k+i) = 1; v(1:n‐k+i‐1) is stored on exit
84in A(m‐k+i,1:n‐k+i‐1), and taua in TAUA(i).
85To form Q explicitly, use LAPACK subroutine ZUNGRQ.
86To use Q to update another matrix, use LAPACK subroutine ZUNMRQ.
87
88The matrix Z is represented as a product of elementary reflectors
89
90   Z = H(1) H(2) . . . H(k), where k = min(p,n).
91
92Each H(i) has the form
93
94   H(i) = I ‐ taub * v * v'
95
96where  taub  is  a complex scalar, and v is a complex vector with
97v(1:i‐1) = 0 and  v(i)  =  1;  v(i+1:p)  is  stored  on  exit  in
98B(i+1:p,i), and taub in TAUB(i).
99To form Z explicitly, use LAPACK subroutine ZUNGQR.
100To use Z to update another matrix, use LAPACK subroutine ZUNMQR.
101
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