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