1 SUBROUTINE ZGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
2 $ LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LDB, LWORK, M, N, P
11 * ..
12 * .. Array Arguments ..
13 COMPLEX*16 A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
14 $ WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGGQRF computes a generalized QR factorization of an N-by-M matrix A
21 * and an N-by-P matrix B:
22 *
23 * A = Q*R, B = Q*T*Z,
24 *
25 * where Q is an N-by-N unitary matrix, Z is a P-by-P unitary matrix,
26 * and R and T assume one of the forms:
27 *
28 * if N >= M, R = ( R11 ) M , or if N < M, R = ( R11 R12 ) N,
29 * ( 0 ) N-M N M-N
30 * M
31 *
32 * where R11 is upper triangular, and
33 *
34 * if N <= P, T = ( 0 T12 ) N, or if N > P, T = ( T11 ) N-P,
35 * P-N N ( T21 ) P
36 * P
37 *
38 * where T12 or T21 is upper triangular.
39 *
40 * In particular, if B is square and nonsingular, the GQR factorization
41 * of A and B implicitly gives the QR factorization of inv(B)*A:
42 *
43 * inv(B)*A = Z**H * (inv(T)*R)
44 *
45 * where inv(B) denotes the inverse of the matrix B, and Z**H denotes the
46 * conjugate transpose of matrix Z.
47 *
48 * Arguments
49 * =========
50 *
51 * N (input) INTEGER
52 * The number of rows of the matrices A and B. N >= 0.
53 *
54 * M (input) INTEGER
55 * The number of columns of the matrix A. M >= 0.
56 *
57 * P (input) INTEGER
58 * The number of columns of the matrix B. P >= 0.
59 *
60 * A (input/output) COMPLEX*16 array, dimension (LDA,M)
61 * On entry, the N-by-M matrix A.
62 * On exit, the elements on and above the diagonal of the array
63 * contain the min(N,M)-by-M upper trapezoidal matrix R (R is
64 * upper triangular if N >= M); the elements below the diagonal,
65 * with the array TAUA, represent the unitary matrix Q as a
66 * product of min(N,M) elementary reflectors (see Further
67 * Details).
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,N).
71 *
72 * TAUA (output) COMPLEX*16 array, dimension (min(N,M))
73 * The scalar factors of the elementary reflectors which
74 * represent the unitary matrix Q (see Further Details).
75 *
76 * B (input/output) COMPLEX*16 array, dimension (LDB,P)
77 * On entry, the N-by-P matrix B.
78 * On exit, if N <= P, the upper triangle of the subarray
79 * B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
80 * if N > P, the elements on and above the (N-P)-th subdiagonal
81 * contain the N-by-P upper trapezoidal matrix T; the remaining
82 * elements, with the array TAUB, represent the unitary
83 * matrix Z as a product of elementary reflectors (see Further
84 * Details).
85 *
86 * LDB (input) INTEGER
87 * The leading dimension of the array B. LDB >= max(1,N).
88 *
89 * TAUB (output) COMPLEX*16 array, dimension (min(N,P))
90 * The scalar factors of the elementary reflectors which
91 * represent the unitary matrix Z (see Further Details).
92 *
93 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
94 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
95 *
96 * LWORK (input) INTEGER
97 * The dimension of the array WORK. LWORK >= max(1,N,M,P).
98 * For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
99 * where NB1 is the optimal blocksize for the QR factorization
100 * of an N-by-M matrix, NB2 is the optimal blocksize for the
101 * RQ factorization of an N-by-P matrix, and NB3 is the optimal
102 * blocksize for a call of ZUNMQR.
103 *
104 * If LWORK = -1, then a workspace query is assumed; the routine
105 * only calculates the optimal size of the WORK array, returns
106 * this value as the first entry of the WORK array, and no error
107 * message related to LWORK is issued by XERBLA.
108 *
109 * INFO (output) INTEGER
110 * = 0: successful exit
111 * < 0: if INFO = -i, the i-th argument had an illegal value.
112 *
113 * Further Details
114 * ===============
115 *
116 * The matrix Q is represented as a product of elementary reflectors
117 *
118 * Q = H(1) H(2) . . . H(k), where k = min(n,m).
119 *
120 * Each H(i) has the form
121 *
122 * H(i) = I - taua * v * v**H
123 *
124 * where taua is a complex scalar, and v is a complex vector with
125 * v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
126 * and taua in TAUA(i).
127 * To form Q explicitly, use LAPACK subroutine ZUNGQR.
128 * To use Q to update another matrix, use LAPACK subroutine ZUNMQR.
129 *
130 * The matrix Z is represented as a product of elementary reflectors
131 *
132 * Z = H(1) H(2) . . . H(k), where k = min(n,p).
133 *
134 * Each H(i) has the form
135 *
136 * H(i) = I - taub * v * v**H
137 *
138 * where taub is a complex scalar, and v is a complex vector with
139 * v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
140 * B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
141 * To form Z explicitly, use LAPACK subroutine ZUNGRQ.
142 * To use Z to update another matrix, use LAPACK subroutine ZUNMRQ.
143 *
144 * =====================================================================
145 *
146 * .. Local Scalars ..
147 LOGICAL LQUERY
148 INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3
149 * ..
150 * .. External Subroutines ..
151 EXTERNAL XERBLA, ZGEQRF, ZGERQF, ZUNMQR
152 * ..
153 * .. External Functions ..
154 INTEGER ILAENV
155 EXTERNAL ILAENV
156 * ..
157 * .. Intrinsic Functions ..
158 INTRINSIC INT, MAX, MIN
159 * ..
160 * .. Executable Statements ..
161 *
162 * Test the input parameters
163 *
164 INFO = 0
165 NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 )
166 NB2 = ILAENV( 1, 'ZGERQF', ' ', N, P, -1, -1 )
167 NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 )
168 NB = MAX( NB1, NB2, NB3 )
169 LWKOPT = MAX( N, M, P )*NB
170 WORK( 1 ) = LWKOPT
171 LQUERY = ( LWORK.EQ.-1 )
172 IF( N.LT.0 ) THEN
173 INFO = -1
174 ELSE IF( M.LT.0 ) THEN
175 INFO = -2
176 ELSE IF( P.LT.0 ) THEN
177 INFO = -3
178 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
179 INFO = -5
180 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
181 INFO = -8
182 ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN
183 INFO = -11
184 END IF
185 IF( INFO.NE.0 ) THEN
186 CALL XERBLA( 'ZGGQRF', -INFO )
187 RETURN
188 ELSE IF( LQUERY ) THEN
189 RETURN
190 END IF
191 *
192 * QR factorization of N-by-M matrix A: A = Q*R
193 *
194 CALL ZGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
195 LOPT = WORK( 1 )
196 *
197 * Update B := Q**H*B.
198 *
199 CALL ZUNMQR( 'Left', 'Conjugate Transpose', N, P, MIN( N, M ), A,
200 $ LDA, TAUA, B, LDB, WORK, LWORK, INFO )
201 LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
202 *
203 * RQ factorization of N-by-P matrix B: B = T*Z.
204 *
205 CALL ZGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
206 WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
207 *
208 RETURN
209 *
210 * End of ZGGQRF
211 *
212 END
2 $ LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LDB, LWORK, M, N, P
11 * ..
12 * .. Array Arguments ..
13 COMPLEX*16 A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
14 $ WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGGQRF computes a generalized QR factorization of an N-by-M matrix A
21 * and an N-by-P matrix B:
22 *
23 * A = Q*R, B = Q*T*Z,
24 *
25 * where Q is an N-by-N unitary matrix, Z is a P-by-P unitary matrix,
26 * and R and T assume one of the forms:
27 *
28 * if N >= M, R = ( R11 ) M , or if N < M, R = ( R11 R12 ) N,
29 * ( 0 ) N-M N M-N
30 * M
31 *
32 * where R11 is upper triangular, and
33 *
34 * if N <= P, T = ( 0 T12 ) N, or if N > P, T = ( T11 ) N-P,
35 * P-N N ( T21 ) P
36 * P
37 *
38 * where T12 or T21 is upper triangular.
39 *
40 * In particular, if B is square and nonsingular, the GQR factorization
41 * of A and B implicitly gives the QR factorization of inv(B)*A:
42 *
43 * inv(B)*A = Z**H * (inv(T)*R)
44 *
45 * where inv(B) denotes the inverse of the matrix B, and Z**H denotes the
46 * conjugate transpose of matrix Z.
47 *
48 * Arguments
49 * =========
50 *
51 * N (input) INTEGER
52 * The number of rows of the matrices A and B. N >= 0.
53 *
54 * M (input) INTEGER
55 * The number of columns of the matrix A. M >= 0.
56 *
57 * P (input) INTEGER
58 * The number of columns of the matrix B. P >= 0.
59 *
60 * A (input/output) COMPLEX*16 array, dimension (LDA,M)
61 * On entry, the N-by-M matrix A.
62 * On exit, the elements on and above the diagonal of the array
63 * contain the min(N,M)-by-M upper trapezoidal matrix R (R is
64 * upper triangular if N >= M); the elements below the diagonal,
65 * with the array TAUA, represent the unitary matrix Q as a
66 * product of min(N,M) elementary reflectors (see Further
67 * Details).
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,N).
71 *
72 * TAUA (output) COMPLEX*16 array, dimension (min(N,M))
73 * The scalar factors of the elementary reflectors which
74 * represent the unitary matrix Q (see Further Details).
75 *
76 * B (input/output) COMPLEX*16 array, dimension (LDB,P)
77 * On entry, the N-by-P matrix B.
78 * On exit, if N <= P, the upper triangle of the subarray
79 * B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
80 * if N > P, the elements on and above the (N-P)-th subdiagonal
81 * contain the N-by-P upper trapezoidal matrix T; the remaining
82 * elements, with the array TAUB, represent the unitary
83 * matrix Z as a product of elementary reflectors (see Further
84 * Details).
85 *
86 * LDB (input) INTEGER
87 * The leading dimension of the array B. LDB >= max(1,N).
88 *
89 * TAUB (output) COMPLEX*16 array, dimension (min(N,P))
90 * The scalar factors of the elementary reflectors which
91 * represent the unitary matrix Z (see Further Details).
92 *
93 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
94 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
95 *
96 * LWORK (input) INTEGER
97 * The dimension of the array WORK. LWORK >= max(1,N,M,P).
98 * For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
99 * where NB1 is the optimal blocksize for the QR factorization
100 * of an N-by-M matrix, NB2 is the optimal blocksize for the
101 * RQ factorization of an N-by-P matrix, and NB3 is the optimal
102 * blocksize for a call of ZUNMQR.
103 *
104 * If LWORK = -1, then a workspace query is assumed; the routine
105 * only calculates the optimal size of the WORK array, returns
106 * this value as the first entry of the WORK array, and no error
107 * message related to LWORK is issued by XERBLA.
108 *
109 * INFO (output) INTEGER
110 * = 0: successful exit
111 * < 0: if INFO = -i, the i-th argument had an illegal value.
112 *
113 * Further Details
114 * ===============
115 *
116 * The matrix Q is represented as a product of elementary reflectors
117 *
118 * Q = H(1) H(2) . . . H(k), where k = min(n,m).
119 *
120 * Each H(i) has the form
121 *
122 * H(i) = I - taua * v * v**H
123 *
124 * where taua is a complex scalar, and v is a complex vector with
125 * v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
126 * and taua in TAUA(i).
127 * To form Q explicitly, use LAPACK subroutine ZUNGQR.
128 * To use Q to update another matrix, use LAPACK subroutine ZUNMQR.
129 *
130 * The matrix Z is represented as a product of elementary reflectors
131 *
132 * Z = H(1) H(2) . . . H(k), where k = min(n,p).
133 *
134 * Each H(i) has the form
135 *
136 * H(i) = I - taub * v * v**H
137 *
138 * where taub is a complex scalar, and v is a complex vector with
139 * v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
140 * B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
141 * To form Z explicitly, use LAPACK subroutine ZUNGRQ.
142 * To use Z to update another matrix, use LAPACK subroutine ZUNMRQ.
143 *
144 * =====================================================================
145 *
146 * .. Local Scalars ..
147 LOGICAL LQUERY
148 INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3
149 * ..
150 * .. External Subroutines ..
151 EXTERNAL XERBLA, ZGEQRF, ZGERQF, ZUNMQR
152 * ..
153 * .. External Functions ..
154 INTEGER ILAENV
155 EXTERNAL ILAENV
156 * ..
157 * .. Intrinsic Functions ..
158 INTRINSIC INT, MAX, MIN
159 * ..
160 * .. Executable Statements ..
161 *
162 * Test the input parameters
163 *
164 INFO = 0
165 NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 )
166 NB2 = ILAENV( 1, 'ZGERQF', ' ', N, P, -1, -1 )
167 NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 )
168 NB = MAX( NB1, NB2, NB3 )
169 LWKOPT = MAX( N, M, P )*NB
170 WORK( 1 ) = LWKOPT
171 LQUERY = ( LWORK.EQ.-1 )
172 IF( N.LT.0 ) THEN
173 INFO = -1
174 ELSE IF( M.LT.0 ) THEN
175 INFO = -2
176 ELSE IF( P.LT.0 ) THEN
177 INFO = -3
178 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
179 INFO = -5
180 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
181 INFO = -8
182 ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN
183 INFO = -11
184 END IF
185 IF( INFO.NE.0 ) THEN
186 CALL XERBLA( 'ZGGQRF', -INFO )
187 RETURN
188 ELSE IF( LQUERY ) THEN
189 RETURN
190 END IF
191 *
192 * QR factorization of N-by-M matrix A: A = Q*R
193 *
194 CALL ZGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
195 LOPT = WORK( 1 )
196 *
197 * Update B := Q**H*B.
198 *
199 CALL ZUNMQR( 'Left', 'Conjugate Transpose', N, P, MIN( N, M ), A,
200 $ LDA, TAUA, B, LDB, WORK, LWORK, INFO )
201 LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
202 *
203 * RQ factorization of N-by-P matrix B: B = T*Z.
204 *
205 CALL ZGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
206 WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
207 *
208 RETURN
209 *
210 * End of ZGGQRF
211 *
212 END