3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SGGRQF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggrqf.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggrqf.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggrqf.f">
21 * SUBROUTINE SGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK,
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDB, LWORK, M, N, P
27 * .. Array Arguments ..
28 * REAL A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
38 *> SGGRQF computes a generalized RQ factorization of an M-by-N matrix A
39 *> and a P-by-N matrix B:
41 *> A = R*Q, B = Z*T*Q,
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:
46 *> if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N,
50 *> where R12 or R21 is upper triangular, and
52 *> if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P,
56 *> where T11 is upper triangular.
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):
61 *> A*inv(B) = (R*inv(T))*Z**T
63 *> where inv(B) denotes the inverse of the matrix B, and Z**T denotes the
64 *> transpose of the matrix Z.
73 *> The number of rows of the matrix A. M >= 0.
79 *> The number of rows of the matrix B. P >= 0.
85 *> The number of columns of the matrices A and B. N >= 0.
90 *> A is REAL 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
104 *> The leading dimension of the array A. LDA >= max(1,M).
109 *> TAUA is REAL array, dimension (min(M,N))
110 *> The scalar factors of the elementary reflectors which
111 *> represent the orthogonal matrix Q (see Further Details).
116 *> B is REAL 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).
128 *> The leading dimension of the array B. LDB >= max(1,P).
133 *> TAUB is REAL array, dimension (min(P,N))
134 *> The scalar factors of the elementary reflectors which
135 *> represent the orthogonal matrix Z (see Further Details).
140 *> WORK is REAL array, dimension (MAX(1,LWORK))
141 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
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 SORMRQ.
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.
163 *> = 0: successful exit
164 *> < 0: if INF0= -i, the i-th argument had an illegal value.
170 *> \author Univ. of Tennessee
171 *> \author Univ. of California Berkeley
172 *> \author Univ. of Colorado Denver
175 *> \date November 2011
177 *> \ingroup realOTHERcomputational
179 *> \par Further Details:
180 * =====================
184 *> The matrix Q is represented as a product of elementary reflectors
186 *> Q = H(1) H(2) . . . H(k), where k = min(m,n).
188 *> Each H(i) has the form
190 *> H(i) = I - taua * v * v**T
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 SORGRQ.
196 *> To use Q to update another matrix, use LAPACK subroutine SORMRQ.
198 *> The matrix Z is represented as a product of elementary reflectors
200 *> Z = H(1) H(2) . . . H(k), where k = min(p,n).
202 *> Each H(i) has the form
204 *> H(i) = I - taub * v * v**T
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 SORGQR.
210 *> To use Z to update another matrix, use LAPACK subroutine SORMQR.
213 * =====================================================================
214 SUBROUTINE SGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK,
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..--
222 * .. Scalar Arguments ..
223 INTEGER INFO, LDA, LDB, LWORK, M, N, P
225 * .. Array Arguments ..
226 REAL A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
230 * =====================================================================
232 * .. Local Scalars ..
234 INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3
236 * .. External Subroutines ..
237 EXTERNAL SGEQRF, SGERQF, SORMRQ, XERBLA
239 * .. External Functions ..
243 * .. Intrinsic Functions ..
244 INTRINSIC INT, MAX, MIN
246 * .. Executable Statements ..
248 * Test the input parameters
251 NB1 = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
252 NB2 = ILAENV( 1, 'SGEQRF', ' ', P, N, -1, -1 )
253 NB3 = ILAENV( 1, 'SORMRQ', ' ', M, N, P, -1 )
254 NB = MAX( NB1, NB2, NB3 )
255 LWKOPT = MAX( N, M, P)*NB
257 LQUERY = ( LWORK.EQ.-1 )
260 ELSE IF( P.LT.0 ) THEN
262 ELSE IF( N.LT.0 ) THEN
264 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
266 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
268 ELSE IF( LWORK.LT.MAX( 1, M, P, N ) .AND. .NOT.LQUERY ) THEN
272 CALL XERBLA( 'SGGRQF', -INFO )
274 ELSE IF( LQUERY ) THEN
278 * RQ factorization of M-by-N matrix A: A = R*Q
280 CALL SGERQF( M, N, A, LDA, TAUA, WORK, LWORK, INFO )
285 CALL SORMRQ( 'Right', 'Transpose', P, N, MIN( M, N ),
286 $ A( MAX( 1, M-N+1 ), 1 ), LDA, TAUA, B, LDB, WORK,
288 LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
290 * QR factorization of P-by-N matrix B: B = Z*T
292 CALL SGEQRF( P, N, B, LDB, TAUB, WORK, LWORK, INFO )
293 WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )