1ZGGQRF ‐ a generalized QR factorization of an N‐by‐M matrix A and
2an N‐by‐P matrix B SUBROUTINE ZGGQRF( N, M, P, 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( * ) ZGGQRF computes a generalized QR factorization of an N‐
7by‐M matrix 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*16  ar‐
39ray,  dimension  (LDA,M) On entry, the N‐by‐M matrix A.  On exit,
40the elements 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*16 array, dimension (min(N,M)) The  scalar  fac‐
47tors of the elementary reflectors which represent the unitary ma‐
48trix Q (see Further Details).  B       (input/output)  COMPLEX*16
49array, dimension (LDB,P) On entry, the N‐by‐P matrix B.  On exit,
50if N <= P, the upper triangle of the subarray B(1:N,P‐N+1:P) con‐
51tains  the  N‐by‐N  upper triangular matrix T; if N > P, the ele‐
52ments on and above the (N‐P)‐th subdiagonal  contain  the  N‐by‐P
53upper  trapezoidal matrix T; the remaining elements, with the ar‐
54ray TAUB, represent the unitary matrix Z as a product of  elemen‐
55tary  reflectors  (see Further Details).  LDB     (input) INTEGER
56The leading dimension of the array  B.  LDB  >=  max(1,N).   TAUB
57(output)  COMPLEX*16  array, dimension (min(N,P)) The scalar fac‐
58tors of the elementary reflectors which represent the unitary ma‐
59trix  Z  (see  Further Details).  WORK    (workspace/output) COM‐
60PLEX*16 array, dimension (MAX(1,LWORK)) On exit,  if  INFO  =  0,
61WORK(1)  returns  the optimal LWORK.  LWORK   (input) INTEGER The
62dimension of the array WORK. LWORK >= max(1,N,M,P).  For  optimum
63performance  LWORK  >=  max(N,M,P)*max(NB1,NB2,NB3), where NB1 is
64the optimal blocksize for the QR factorization of an  N‐by‐M  ma‐
65trix, NB2 is the optimal blocksize for the RQ factorization of an
66N‐by‐P matrix, and NB3 is the optimal blocksize  for  a  call  of
67ZUNMQR.
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 INFO = ‐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(n,m).
78
79Each H(i) has the form
80
81   H(i) = I ‐ taua * v * v'
82
83where taua is a complex scalar, and v is a  complex  vector  with
84v(1:i‐1)  =  0  and  v(i)  =  1;  v(i+1:n)  is  stored on exit in
85A(i+1:n,i), and taua in TAUA(i).
86To form Q explicitly, use LAPACK subroutine ZUNGQR.
87To use Q to update another matrix, use LAPACK subroutine ZUNMQR.
88
89The matrix Z is represented as a product of elementary reflectors
90
91   Z = H(1) H(2) . . . H(k), where k = min(n,p).
92
93Each H(i) has the form
94
95   H(i) = I ‐ taub * v * v'
96
97where taub is a complex scalar, and v is a  complex  vector  with
98v(p‐k+i+1:p) = 0 and v(p‐k+i) = 1; v(1:p‐k+i‐1) is stored on exit
99in B(n‐k+i,1:p‐k+i‐1), and taub in TAUB(i).
100To form Z explicitly, use LAPACK subroutine ZUNGRQ.
101To use Z to update another matrix, use LAPACK subroutine ZUNMRQ.
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