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