![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGGRQF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DGGRQF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggrqf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggrqf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggrqf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, 00022 * LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, LWORK, M, N, P 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ), 00029 * $ WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DGGRQF computes a generalized RQ factorization of an M-by-N matrix A 00039 *> and a P-by-N matrix B: 00040 *> 00041 *> A = R*Q, B = Z*T*Q, 00042 *> 00043 *> where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal 00044 *> matrix, and R and T assume one of the forms: 00045 *> 00046 *> if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N, 00047 *> N-M M ( R21 ) N 00048 *> N 00049 *> 00050 *> where R12 or R21 is upper triangular, and 00051 *> 00052 *> if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P, 00053 *> ( 0 ) P-N P N-P 00054 *> N 00055 *> 00056 *> where T11 is upper triangular. 00057 *> 00058 *> In particular, if B is square and nonsingular, the GRQ factorization 00059 *> of A and B implicitly gives the RQ factorization of A*inv(B): 00060 *> 00061 *> A*inv(B) = (R*inv(T))*Z**T 00062 *> 00063 *> where inv(B) denotes the inverse of the matrix B, and Z**T denotes the 00064 *> transpose of the matrix Z. 00065 *> \endverbatim 00066 * 00067 * Arguments: 00068 * ========== 00069 * 00070 *> \param[in] M 00071 *> \verbatim 00072 *> M is INTEGER 00073 *> The number of rows of the matrix A. M >= 0. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] P 00077 *> \verbatim 00078 *> P is INTEGER 00079 *> The number of rows of the matrix B. P >= 0. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] N 00083 *> \verbatim 00084 *> N is INTEGER 00085 *> The number of columns of the matrices A and B. N >= 0. 00086 *> \endverbatim 00087 *> 00088 *> \param[in,out] A 00089 *> \verbatim 00090 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00091 *> On entry, the M-by-N matrix A. 00092 *> On exit, if M <= N, the upper triangle of the subarray 00093 *> A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R; 00094 *> if M > N, the elements on and above the (M-N)-th subdiagonal 00095 *> contain the M-by-N upper trapezoidal matrix R; the remaining 00096 *> elements, with the array TAUA, represent the orthogonal 00097 *> matrix Q as a product of elementary reflectors (see Further 00098 *> Details). 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDA 00102 *> \verbatim 00103 *> LDA is INTEGER 00104 *> The leading dimension of the array A. LDA >= max(1,M). 00105 *> \endverbatim 00106 *> 00107 *> \param[out] TAUA 00108 *> \verbatim 00109 *> TAUA is DOUBLE PRECISION array, dimension (min(M,N)) 00110 *> The scalar factors of the elementary reflectors which 00111 *> represent the orthogonal matrix Q (see Further Details). 00112 *> \endverbatim 00113 *> 00114 *> \param[in,out] B 00115 *> \verbatim 00116 *> B is DOUBLE PRECISION array, dimension (LDB,N) 00117 *> On entry, the P-by-N matrix B. 00118 *> On exit, the elements on and above the diagonal of the array 00119 *> contain the min(P,N)-by-N upper trapezoidal matrix T (T is 00120 *> upper triangular if P >= N); the elements below the diagonal, 00121 *> with the array TAUB, represent the orthogonal matrix Z as a 00122 *> product of elementary reflectors (see Further Details). 00123 *> \endverbatim 00124 *> 00125 *> \param[in] LDB 00126 *> \verbatim 00127 *> LDB is INTEGER 00128 *> The leading dimension of the array B. LDB >= max(1,P). 00129 *> \endverbatim 00130 *> 00131 *> \param[out] TAUB 00132 *> \verbatim 00133 *> TAUB is DOUBLE PRECISION array, dimension (min(P,N)) 00134 *> The scalar factors of the elementary reflectors which 00135 *> represent the orthogonal matrix Z (see Further Details). 00136 *> \endverbatim 00137 *> 00138 *> \param[out] WORK 00139 *> \verbatim 00140 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00141 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00142 *> \endverbatim 00143 *> 00144 *> \param[in] LWORK 00145 *> \verbatim 00146 *> LWORK is INTEGER 00147 *> The dimension of the array WORK. LWORK >= max(1,N,M,P). 00148 *> For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3), 00149 *> where NB1 is the optimal blocksize for the RQ factorization 00150 *> of an M-by-N matrix, NB2 is the optimal blocksize for the 00151 *> QR factorization of a P-by-N matrix, and NB3 is the optimal 00152 *> blocksize for a call of DORMRQ. 00153 *> 00154 *> If LWORK = -1, then a workspace query is assumed; the routine 00155 *> only calculates the optimal size of the WORK array, returns 00156 *> this value as the first entry of the WORK array, and no error 00157 *> message related to LWORK is issued by XERBLA. 00158 *> \endverbatim 00159 *> 00160 *> \param[out] INFO 00161 *> \verbatim 00162 *> INFO is INTEGER 00163 *> = 0: successful exit 00164 *> < 0: if INF0= -i, the i-th argument had an illegal value. 00165 *> \endverbatim 00166 * 00167 * Authors: 00168 * ======== 00169 * 00170 *> \author Univ. of Tennessee 00171 *> \author Univ. of California Berkeley 00172 *> \author Univ. of Colorado Denver 00173 *> \author NAG Ltd. 00174 * 00175 *> \date November 2011 00176 * 00177 *> \ingroup doubleOTHERcomputational 00178 * 00179 *> \par Further Details: 00180 * ===================== 00181 *> 00182 *> \verbatim 00183 *> 00184 *> The matrix Q is represented as a product of elementary reflectors 00185 *> 00186 *> Q = H(1) H(2) . . . H(k), where k = min(m,n). 00187 *> 00188 *> Each H(i) has the form 00189 *> 00190 *> H(i) = I - taua * v * v**T 00191 *> 00192 *> where taua is a real scalar, and v is a real vector with 00193 *> v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in 00194 *> A(m-k+i,1:n-k+i-1), and taua in TAUA(i). 00195 *> To form Q explicitly, use LAPACK subroutine DORGRQ. 00196 *> To use Q to update another matrix, use LAPACK subroutine DORMRQ. 00197 *> 00198 *> The matrix Z is represented as a product of elementary reflectors 00199 *> 00200 *> Z = H(1) H(2) . . . H(k), where k = min(p,n). 00201 *> 00202 *> Each H(i) has the form 00203 *> 00204 *> H(i) = I - taub * v * v**T 00205 *> 00206 *> where taub is a real scalar, and v is a real vector with 00207 *> v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i), 00208 *> and taub in TAUB(i). 00209 *> To form Z explicitly, use LAPACK subroutine DORGQR. 00210 *> To use Z to update another matrix, use LAPACK subroutine DORMQR. 00211 *> \endverbatim 00212 *> 00213 * ===================================================================== 00214 SUBROUTINE DGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, 00215 $ LWORK, INFO ) 00216 * 00217 * -- LAPACK computational routine (version 3.4.0) -- 00218 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00220 * November 2011 00221 * 00222 * .. Scalar Arguments .. 00223 INTEGER INFO, LDA, LDB, LWORK, M, N, P 00224 * .. 00225 * .. Array Arguments .. 00226 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ), 00227 $ WORK( * ) 00228 * .. 00229 * 00230 * ===================================================================== 00231 * 00232 * .. Local Scalars .. 00233 LOGICAL LQUERY 00234 INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3 00235 * .. 00236 * .. External Subroutines .. 00237 EXTERNAL DGEQRF, DGERQF, DORMRQ, XERBLA 00238 * .. 00239 * .. External Functions .. 00240 INTEGER ILAENV 00241 EXTERNAL ILAENV 00242 * .. 00243 * .. Intrinsic Functions .. 00244 INTRINSIC INT, MAX, MIN 00245 * .. 00246 * .. Executable Statements .. 00247 * 00248 * Test the input parameters 00249 * 00250 INFO = 0 00251 NB1 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 ) 00252 NB2 = ILAENV( 1, 'DGEQRF', ' ', P, N, -1, -1 ) 00253 NB3 = ILAENV( 1, 'DORMRQ', ' ', M, N, P, -1 ) 00254 NB = MAX( NB1, NB2, NB3 ) 00255 LWKOPT = MAX( N, M, P )*NB 00256 WORK( 1 ) = LWKOPT 00257 LQUERY = ( LWORK.EQ.-1 ) 00258 IF( M.LT.0 ) THEN 00259 INFO = -1 00260 ELSE IF( P.LT.0 ) THEN 00261 INFO = -2 00262 ELSE IF( N.LT.0 ) THEN 00263 INFO = -3 00264 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00265 INFO = -5 00266 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00267 INFO = -8 00268 ELSE IF( LWORK.LT.MAX( 1, M, P, N ) .AND. .NOT.LQUERY ) THEN 00269 INFO = -11 00270 END IF 00271 IF( INFO.NE.0 ) THEN 00272 CALL XERBLA( 'DGGRQF', -INFO ) 00273 RETURN 00274 ELSE IF( LQUERY ) THEN 00275 RETURN 00276 END IF 00277 * 00278 * RQ factorization of M-by-N matrix A: A = R*Q 00279 * 00280 CALL DGERQF( M, N, A, LDA, TAUA, WORK, LWORK, INFO ) 00281 LOPT = WORK( 1 ) 00282 * 00283 * Update B := B*Q**T 00284 * 00285 CALL DORMRQ( 'Right', 'Transpose', P, N, MIN( M, N ), 00286 $ A( MAX( 1, M-N+1 ), 1 ), LDA, TAUA, B, LDB, WORK, 00287 $ LWORK, INFO ) 00288 LOPT = MAX( LOPT, INT( WORK( 1 ) ) ) 00289 * 00290 * QR factorization of P-by-N matrix B: B = Z*T 00291 * 00292 CALL DGEQRF( P, N, B, LDB, TAUB, WORK, LWORK, INFO ) 00293 WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) ) 00294 * 00295 RETURN 00296 * 00297 * End of DGGRQF 00298 * 00299 END