LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dggrqf.f
Go to the documentation of this file.
1 *> \brief \b DGGRQF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DGGRQF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggrqf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggrqf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggrqf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK,
22 * LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDB, LWORK, M, N, P
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
29 * $ WORK( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DGGRQF computes a generalized RQ factorization of an M-by-N matrix A
39 *> and a P-by-N matrix B:
40 *>
41 *> A = R*Q, B = Z*T*Q,
42 *>
43 *> where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
44 *> matrix, and R and T assume one of the forms:
45 *>
46 *> if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N,
47 *> N-M M ( R21 ) N
48 *> N
49 *>
50 *> where R12 or R21 is upper triangular, and
51 *>
52 *> if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P,
53 *> ( 0 ) P-N P N-P
54 *> N
55 *>
56 *> where T11 is upper triangular.
57 *>
58 *> In particular, if B is square and nonsingular, the GRQ factorization
59 *> of A and B implicitly gives the RQ factorization of A*inv(B):
60 *>
61 *> A*inv(B) = (R*inv(T))*Z**T
62 *>
63 *> where inv(B) denotes the inverse of the matrix B, and Z**T denotes the
64 *> transpose of the matrix Z.
65 *> \endverbatim
66 *
67 * Arguments:
68 * ==========
69 *
70 *> \param[in] M
71 *> \verbatim
72 *> M is INTEGER
73 *> The number of rows of the matrix A. M >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in] P
77 *> \verbatim
78 *> P is INTEGER
79 *> The number of rows of the matrix B. P >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in] N
83 *> \verbatim
84 *> N is INTEGER
85 *> The number of columns of the matrices A and B. N >= 0.
86 *> \endverbatim
87 *>
88 *> \param[in,out] A
89 *> \verbatim
90 *> A is DOUBLE PRECISION array, dimension (LDA,N)
91 *> On entry, the M-by-N matrix A.
92 *> On exit, if M <= N, the upper triangle of the subarray
93 *> A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R;
94 *> if M > N, the elements on and above the (M-N)-th subdiagonal
95 *> contain the M-by-N upper trapezoidal matrix R; the remaining
96 *> elements, with the array TAUA, represent the orthogonal
97 *> matrix Q as a product of elementary reflectors (see Further
98 *> Details).
99 *> \endverbatim
100 *>
101 *> \param[in] LDA
102 *> \verbatim
103 *> LDA is INTEGER
104 *> The leading dimension of the array A. LDA >= max(1,M).
105 *> \endverbatim
106 *>
107 *> \param[out] TAUA
108 *> \verbatim
109 *> TAUA is DOUBLE PRECISION array, dimension (min(M,N))
110 *> The scalar factors of the elementary reflectors which
111 *> represent the orthogonal matrix Q (see Further Details).
112 *> \endverbatim
113 *>
114 *> \param[in,out] B
115 *> \verbatim
116 *> B is DOUBLE PRECISION array, dimension (LDB,N)
117 *> On entry, the P-by-N matrix B.
118 *> On exit, the elements on and above the diagonal of the array
119 *> contain the min(P,N)-by-N upper trapezoidal matrix T (T is
120 *> upper triangular if P >= N); the elements below the diagonal,
121 *> with the array TAUB, represent the orthogonal matrix Z as a
122 *> product of elementary reflectors (see Further Details).
123 *> \endverbatim
124 *>
125 *> \param[in] LDB
126 *> \verbatim
127 *> LDB is INTEGER
128 *> The leading dimension of the array B. LDB >= max(1,P).
129 *> \endverbatim
130 *>
131 *> \param[out] TAUB
132 *> \verbatim
133 *> TAUB is DOUBLE PRECISION array, dimension (min(P,N))
134 *> The scalar factors of the elementary reflectors which
135 *> represent the orthogonal matrix Z (see Further Details).
136 *> \endverbatim
137 *>
138 *> \param[out] WORK
139 *> \verbatim
140 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
141 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
142 *> \endverbatim
143 *>
144 *> \param[in] LWORK
145 *> \verbatim
146 *> LWORK is INTEGER
147 *> The dimension of the array WORK. LWORK >= max(1,N,M,P).
148 *> For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
149 *> where NB1 is the optimal blocksize for the RQ factorization
150 *> of an M-by-N matrix, NB2 is the optimal blocksize for the
151 *> QR factorization of a P-by-N matrix, and NB3 is the optimal
152 *> blocksize for a call of DORMRQ.
153 *>
154 *> If LWORK = -1, then a workspace query is assumed; the routine
155 *> only calculates the optimal size of the WORK array, returns
156 *> this value as the first entry of the WORK array, and no error
157 *> message related to LWORK is issued by XERBLA.
158 *> \endverbatim
159 *>
160 *> \param[out] INFO
161 *> \verbatim
162 *> INFO is INTEGER
163 *> = 0: successful exit
164 *> < 0: if INF0= -i, the i-th argument had an illegal value.
165 *> \endverbatim
166 *
167 * Authors:
168 * ========
169 *
170 *> \author Univ. of Tennessee
171 *> \author Univ. of California Berkeley
172 *> \author Univ. of Colorado Denver
173 *> \author NAG Ltd.
174 *
175 *> \date November 2011
176 *
177 *> \ingroup doubleOTHERcomputational
178 *
179 *> \par Further Details:
180 * =====================
181 *>
182 *> \verbatim
183 *>
184 *> The matrix Q is represented as a product of elementary reflectors
185 *>
186 *> Q = H(1) H(2) . . . H(k), where k = min(m,n).
187 *>
188 *> Each H(i) has the form
189 *>
190 *> H(i) = I - taua * v * v**T
191 *>
192 *> where taua is a real scalar, and v is a real vector with
193 *> v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
194 *> A(m-k+i,1:n-k+i-1), and taua in TAUA(i).
195 *> To form Q explicitly, use LAPACK subroutine DORGRQ.
196 *> To use Q to update another matrix, use LAPACK subroutine DORMRQ.
197 *>
198 *> The matrix Z is represented as a product of elementary reflectors
199 *>
200 *> Z = H(1) H(2) . . . H(k), where k = min(p,n).
201 *>
202 *> Each H(i) has the form
203 *>
204 *> H(i) = I - taub * v * v**T
205 *>
206 *> where taub is a real scalar, and v is a real vector with
207 *> v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i),
208 *> and taub in TAUB(i).
209 *> To form Z explicitly, use LAPACK subroutine DORGQR.
210 *> To use Z to update another matrix, use LAPACK subroutine DORMQR.
211 *> \endverbatim
212 *>
213 * =====================================================================
214  SUBROUTINE dggrqf( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK,
215  $ lwork, info )
216 *
217 * -- LAPACK computational routine (version 3.4.0) --
218 * -- LAPACK is a software package provided by Univ. of Tennessee, --
219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220 * November 2011
221 *
222 * .. Scalar Arguments ..
223  INTEGER info, lda, ldb, lwork, m, n, p
224 * ..
225 * .. Array Arguments ..
226  DOUBLE PRECISION a( lda, * ), b( ldb, * ), taua( * ), taub( * ),
227  $ work( * )
228 * ..
229 *
230 * =====================================================================
231 *
232 * .. Local Scalars ..
233  LOGICAL lquery
234  INTEGER lopt, lwkopt, nb, nb1, nb2, nb3
235 * ..
236 * .. External Subroutines ..
237  EXTERNAL dgeqrf, dgerqf, dormrq, xerbla
238 * ..
239 * .. External Functions ..
240  INTEGER ilaenv
241  EXTERNAL ilaenv
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC int, max, min
245 * ..
246 * .. Executable Statements ..
247 *
248 * Test the input parameters
249 *
250  info = 0
251  nb1 = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
252  nb2 = ilaenv( 1, 'DGEQRF', ' ', p, n, -1, -1 )
253  nb3 = ilaenv( 1, 'DORMRQ', ' ', m, n, p, -1 )
254  nb = max( nb1, nb2, nb3 )
255  lwkopt = max( n, m, p )*nb
256  work( 1 ) = lwkopt
257  lquery = ( lwork.EQ.-1 )
258  IF( m.LT.0 ) THEN
259  info = -1
260  ELSE IF( p.LT.0 ) THEN
261  info = -2
262  ELSE IF( n.LT.0 ) THEN
263  info = -3
264  ELSE IF( lda.LT.max( 1, m ) ) THEN
265  info = -5
266  ELSE IF( ldb.LT.max( 1, p ) ) THEN
267  info = -8
268  ELSE IF( lwork.LT.max( 1, m, p, n ) .AND. .NOT.lquery ) THEN
269  info = -11
270  END IF
271  IF( info.NE.0 ) THEN
272  CALL xerbla( 'DGGRQF', -info )
273  return
274  ELSE IF( lquery ) THEN
275  return
276  END IF
277 *
278 * RQ factorization of M-by-N matrix A: A = R*Q
279 *
280  CALL dgerqf( m, n, a, lda, taua, work, lwork, info )
281  lopt = work( 1 )
282 *
283 * Update B := B*Q**T
284 *
285  CALL dormrq( 'Right', 'Transpose', p, n, min( m, n ),
286  $ a( max( 1, m-n+1 ), 1 ), lda, taua, b, ldb, work,
287  $ lwork, info )
288  lopt = max( lopt, int( work( 1 ) ) )
289 *
290 * QR factorization of P-by-N matrix B: B = Z*T
291 *
292  CALL dgeqrf( p, n, b, ldb, taub, work, lwork, info )
293  work( 1 ) = max( lopt, int( work( 1 ) ) )
294 *
295  return
296 *
297 * End of DGGRQF
298 *
299  END