1CGGQRF ‐ a generalized QR factorization of an N‐by‐M matrix A and
2an N‐by‐P matrix B SUBROUTINE CGGQRF( N, M, P, A, LDA,  TAUA,  B,
3LDB, TAUB, WORK, LWORK, INFO )
4    INTEGER INFO, LDA, LDB, LWORK, M, N, P
5    COMPLEX A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ), WORK(
6* ) CGGQRF computes a generalized QR factorization of  an  N‐by‐M
7matrix A and an N‐by‐P matrix B:
8
9            A = Q*R,        B = Q*T*Z,
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 N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N,
15                (  0  ) N‐M                         N   M‐N
16                   M
17
18where R11 is upper triangular, and
19
20if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N‐P,
21                 P‐N  N                           ( T21 ) P
22                                                     P
23
24where T12 or T21 is upper triangular.
25
26In particular, if B is square and nonsingular, the GQR factoriza‐
27tion  of  A  and  B  implicitly  gives  the  QR  factorization of
28inv(B)*A:
29
30             inv(B)*A = Z'*(inv(T)*R)
31
32where inv(B) denotes the inverse of the matrix B, and Z'  denotes
33the conjugate transpose of matrix Z.
34
35N        (input) INTEGER The number of rows of the matrices A and
36B. N >= 0.  M       (input) INTEGER The number of columns of  the
37matrix A.  M >= 0.  P       (input) INTEGER The number of columns
38of the matrix B.  P >= 0.  A       (input/output) COMPLEX  array,
39dimension  (LDA,M)  On  entry, the N‐by‐M matrix A.  On exit, the
40elements on and above the  diagonal  of  the  array  contain  the
41min(N,M)‐by‐M  upper  trapezoidal matrix R (R is upper triangular
42if N >= M); the elements below the diagonal, with the array TAUA,
43represent  the  unitary matrix Q as a product of min(N,M) elemen‐
44tary reflectors (see Further Details).  LDA      (input)  INTEGER
45The  leading  dimension  of  the  array A. LDA >= max(1,N).  TAUA
46(output) COMPLEX array, dimension (min(N,M)) The  scalar  factors
47of the elementary reflectors which represent the unitary matrix Q
48(see Further Details).  B       (input/output) COMPLEX array, di‐
49mension  (LDB,P) On entry, the N‐by‐P matrix B.  On exit, if N <=
50P, the upper triangle of the subarray B(1:N,P‐N+1:P) contains the
51N‐by‐N  upper  triangular matrix T; if N > P, the elements on and
52above the (N‐P)‐th subdiagonal contain the  N‐by‐P  upper  trape‐
53zoidal  matrix  T;  the  remaining elements, 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,N).  TAUB    (output) COM‐
57PLEX  array,  dimension (min(N,P)) The scalar factors of the ele‐
58mentary reflectors which represent the unitary matrix Z (see Fur‐
59ther  Details).  WORK    (workspace/output) COMPLEX array, dimen‐
60sion (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the op‐
61timal  LWORK.  LWORK   (input) INTEGER The dimension of the array
62WORK. LWORK >= max(1,N,M,P).  For optimum  performance  LWORK  >=
63max(N,M,P)*max(NB1,NB2,NB3),  where  NB1 is the optimal blocksize
64for the QR factorization of an N‐by‐M matrix, NB2 is the  optimal
65blocksize  for  the RQ factorization of an N‐by‐P matrix, and NB3
66is the optimal blocksize for a call of CUNMQR.
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(n,m).
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(1:i‐1) = 0 and  v(i)  =  1;  v(i+1:n)  is  stored  on  exit  in
84A(i+1:n,i), and taua in TAUA(i).
85To form Q explicitly, use LAPACK subroutine CUNGQR.
86To use Q to update another matrix, use LAPACK subroutine CUNMQR.
87
88The matrix Z is represented as a product of elementary reflectors
89
90   Z = H(1) H(2) . . . H(k), where k = min(n,p).
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(p‐k+i+1:p) = 0 and v(p‐k+i) = 1; v(1:p‐k+i‐1) is stored on exit
98in B(n‐k+i,1:p‐k+i‐1), and taub in TAUB(i).
99To form Z explicitly, use LAPACK subroutine CUNGRQ.
100To use Z to update another matrix, use LAPACK subroutine CUNMRQ.
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