![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 C> \brief \b SGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm. 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE SGEQRF ( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER INFO, LDA, LWORK, M, N 00015 * .. 00016 * .. Array Arguments .. 00017 * REAL A( LDA, * ), TAU( * ), WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 C>\details \b Purpose: 00024 C>\verbatim 00025 C> 00026 C> SGEQRF computes a QR factorization of a real M-by-N matrix A: 00027 C> A = Q * R. 00028 C> 00029 C> This is the left-looking Level 3 BLAS version of the algorithm. 00030 C> 00031 C>\endverbatim 00032 * 00033 * Arguments: 00034 * ========== 00035 * 00036 C> \param[in] M 00037 C> \verbatim 00038 C> M is INTEGER 00039 C> The number of rows of the matrix A. M >= 0. 00040 C> \endverbatim 00041 C> 00042 C> \param[in] N 00043 C> \verbatim 00044 C> N is INTEGER 00045 C> The number of columns of the matrix A. N >= 0. 00046 C> \endverbatim 00047 C> 00048 C> \param[in,out] A 00049 C> \verbatim 00050 C> A is REAL array, dimension (LDA,N) 00051 C> On entry, the M-by-N matrix A. 00052 C> On exit, the elements on and above the diagonal of the array 00053 C> contain the min(M,N)-by-N upper trapezoidal matrix R (R is 00054 C> upper triangular if m >= n); the elements below the diagonal, 00055 C> with the array TAU, represent the orthogonal matrix Q as a 00056 C> product of min(m,n) elementary reflectors (see Further 00057 C> Details). 00058 C> \endverbatim 00059 C> 00060 C> \param[in] LDA 00061 C> \verbatim 00062 C> LDA is INTEGER 00063 C> The leading dimension of the array A. LDA >= max(1,M). 00064 C> \endverbatim 00065 C> 00066 C> \param[out] TAU 00067 C> \verbatim 00068 C> TAU is REAL array, dimension (min(M,N)) 00069 C> The scalar factors of the elementary reflectors (see Further 00070 C> Details). 00071 C> \endverbatim 00072 C> 00073 C> \param[out] WORK 00074 C> \verbatim 00075 C> WORK is REAL array, dimension (MAX(1,LWORK)) 00076 C> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00077 C> \endverbatim 00078 C> 00079 C> \param[in] LWORK 00080 C> \verbatim 00081 C> LWORK is INTEGER 00082 C> \endverbatim 00083 C> \verbatim 00084 C> The dimension of the array WORK. The dimension can be divided into three parts. 00085 C> \endverbatim 00086 C> \verbatim 00087 C> 1) The part for the triangular factor T. If the very last T is not bigger 00088 C> than any of the rest, then this part is NB x ceiling(K/NB), otherwise, 00089 C> NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T 00090 C> \endverbatim 00091 C> \verbatim 00092 C> 2) The part for the very last T when T is bigger than any of the rest T. 00093 C> The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB, 00094 C> where K = min(M,N), NX is calculated by 00095 C> NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) ) 00096 C> \endverbatim 00097 C> \verbatim 00098 C> 3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB) 00099 C> \endverbatim 00100 C> \verbatim 00101 C> So LWORK = part1 + part2 + part3 00102 C> \endverbatim 00103 C> \verbatim 00104 C> If LWORK = -1, then a workspace query is assumed; the routine 00105 C> only calculates the optimal size of the WORK array, returns 00106 C> this value as the first entry of the WORK array, and no error 00107 C> message related to LWORK is issued by XERBLA. 00108 C> \endverbatim 00109 C> 00110 C> \param[out] INFO 00111 C> \verbatim 00112 C> INFO is INTEGER 00113 C> = 0: successful exit 00114 C> < 0: if INFO = -i, the i-th argument had an illegal value 00115 C> \endverbatim 00116 C> 00117 * 00118 * Authors: 00119 * ======== 00120 * 00121 C> \author Univ. of Tennessee 00122 C> \author Univ. of California Berkeley 00123 C> \author Univ. of Colorado Denver 00124 C> \author NAG Ltd. 00125 * 00126 C> \date November 2011 00127 * 00128 C> \ingroup variantsGEcomputational 00129 * 00130 * Further Details 00131 * =============== 00132 C>\details \b Further \b Details 00133 C> \verbatim 00134 C> 00135 C> The matrix Q is represented as a product of elementary reflectors 00136 C> 00137 C> Q = H(1) H(2) . . . H(k), where k = min(m,n). 00138 C> 00139 C> Each H(i) has the form 00140 C> 00141 C> H(i) = I - tau * v * v' 00142 C> 00143 C> where tau is a real scalar, and v is a real vector with 00144 C> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), 00145 C> and tau in TAU(i). 00146 C> 00147 C> \endverbatim 00148 C> 00149 * ===================================================================== 00150 SUBROUTINE SGEQRF ( M, N, A, LDA, TAU, WORK, LWORK, INFO ) 00151 * 00152 * -- LAPACK computational routine (version 3.1) -- 00153 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00155 * November 2011 00156 * 00157 * .. Scalar Arguments .. 00158 INTEGER INFO, LDA, LWORK, M, N 00159 * .. 00160 * .. Array Arguments .. 00161 REAL A( LDA, * ), TAU( * ), WORK( * ) 00162 * .. 00163 * 00164 * ===================================================================== 00165 * 00166 * .. Local Scalars .. 00167 LOGICAL LQUERY 00168 INTEGER I, IB, IINFO, IWS, J, K, LWKOPT, NB, 00169 $ NBMIN, NX, LBWORK, NT, LLWORK 00170 * .. 00171 * .. External Subroutines .. 00172 EXTERNAL SGEQR2, SLARFB, SLARFT, XERBLA 00173 * .. 00174 * .. Intrinsic Functions .. 00175 INTRINSIC MAX, MIN 00176 * .. 00177 * .. External Functions .. 00178 INTEGER ILAENV 00179 REAL SCEIL 00180 EXTERNAL ILAENV, SCEIL 00181 * .. 00182 * .. Executable Statements .. 00183 00184 INFO = 0 00185 NBMIN = 2 00186 NX = 0 00187 IWS = N 00188 K = MIN( M, N ) 00189 NB = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 ) 00190 00191 IF( NB.GT.1 .AND. NB.LT.K ) THEN 00192 * 00193 * Determine when to cross over from blocked to unblocked code. 00194 * 00195 NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) ) 00196 END IF 00197 * 00198 * Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.: 00199 * 00200 * NB=3 2NB=6 K=10 00201 * | | | 00202 * 1--2--3--4--5--6--7--8--9--10 00203 * | \________/ 00204 * K-NX=5 NT=4 00205 * 00206 * So here 4 x 4 is the last T stored in the workspace 00207 * 00208 NT = K-SCEIL(REAL(K-NX)/REAL(NB))*NB 00209 00210 * 00211 * optimal workspace = space for dlarfb + space for normal T's + space for the last T 00212 * 00213 LLWORK = MAX (MAX((N-M)*K, (N-M)*NB), MAX(K*NB, NB*NB)) 00214 LLWORK = SCEIL(REAL(LLWORK)/REAL(NB)) 00215 00216 IF ( NT.GT.NB ) THEN 00217 00218 LBWORK = K-NT 00219 * 00220 * Optimal workspace for dlarfb = MAX(1,N)*NT 00221 * 00222 LWKOPT = (LBWORK+LLWORK)*NB 00223 WORK( 1 ) = (LWKOPT+NT*NT) 00224 00225 ELSE 00226 00227 LBWORK = SCEIL(REAL(K)/REAL(NB))*NB 00228 LWKOPT = (LBWORK+LLWORK-NB)*NB 00229 WORK( 1 ) = LWKOPT 00230 00231 END IF 00232 00233 * 00234 * Test the input arguments 00235 * 00236 LQUERY = ( LWORK.EQ.-1 ) 00237 IF( M.LT.0 ) THEN 00238 INFO = -1 00239 ELSE IF( N.LT.0 ) THEN 00240 INFO = -2 00241 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00242 INFO = -4 00243 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00244 INFO = -7 00245 END IF 00246 IF( INFO.NE.0 ) THEN 00247 CALL XERBLA( 'SGEQRF', -INFO ) 00248 RETURN 00249 ELSE IF( LQUERY ) THEN 00250 RETURN 00251 END IF 00252 * 00253 * Quick return if possible 00254 * 00255 IF( K.EQ.0 ) THEN 00256 WORK( 1 ) = 1 00257 RETURN 00258 END IF 00259 * 00260 IF( NB.GT.1 .AND. NB.LT.K ) THEN 00261 00262 IF( NX.LT.K ) THEN 00263 * 00264 * Determine if workspace is large enough for blocked code. 00265 * 00266 IF ( NT.LE.NB ) THEN 00267 IWS = (LBWORK+LLWORK-NB)*NB 00268 ELSE 00269 IWS = (LBWORK+LLWORK)*NB+NT*NT 00270 END IF 00271 00272 IF( LWORK.LT.IWS ) THEN 00273 * 00274 * Not enough workspace to use optimal NB: reduce NB and 00275 * determine the minimum value of NB. 00276 * 00277 IF ( NT.LE.NB ) THEN 00278 NB = LWORK / (LLWORK+(LBWORK-NB)) 00279 ELSE 00280 NB = (LWORK-NT*NT)/(LBWORK+LLWORK) 00281 END IF 00282 00283 NBMIN = MAX( 2, ILAENV( 2, 'SGEQRF', ' ', M, N, -1, 00284 $ -1 ) ) 00285 END IF 00286 END IF 00287 END IF 00288 * 00289 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 00290 * 00291 * Use blocked code initially 00292 * 00293 DO 10 I = 1, K - NX, NB 00294 IB = MIN( K-I+1, NB ) 00295 * 00296 * Update the current column using old T's 00297 * 00298 DO 20 J = 1, I - NB, NB 00299 * 00300 * Apply H' to A(J:M,I:I+IB-1) from the left 00301 * 00302 CALL SLARFB( 'Left', 'Transpose', 'Forward', 00303 $ 'Columnwise', M-J+1, IB, NB, 00304 $ A( J, J ), LDA, WORK(J), LBWORK, 00305 $ A( J, I ), LDA, WORK(LBWORK*NB+NT*NT+1), 00306 $ IB) 00307 00308 20 CONTINUE 00309 * 00310 * Compute the QR factorization of the current block 00311 * A(I:M,I:I+IB-1) 00312 * 00313 CALL SGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ), 00314 $ WORK(LBWORK*NB+NT*NT+1), IINFO ) 00315 00316 IF( I+IB.LE.N ) THEN 00317 * 00318 * Form the triangular factor of the block reflector 00319 * H = H(i) H(i+1) . . . H(i+ib-1) 00320 * 00321 CALL SLARFT( 'Forward', 'Columnwise', M-I+1, IB, 00322 $ A( I, I ), LDA, TAU( I ), 00323 $ WORK(I), LBWORK ) 00324 * 00325 END IF 00326 10 CONTINUE 00327 ELSE 00328 I = 1 00329 END IF 00330 * 00331 * Use unblocked code to factor the last or only block. 00332 * 00333 IF( I.LE.K ) THEN 00334 00335 IF ( I .NE. 1 ) THEN 00336 00337 DO 30 J = 1, I - NB, NB 00338 * 00339 * Apply H' to A(J:M,I:K) from the left 00340 * 00341 CALL SLARFB( 'Left', 'Transpose', 'Forward', 00342 $ 'Columnwise', M-J+1, K-I+1, NB, 00343 $ A( J, J ), LDA, WORK(J), LBWORK, 00344 $ A( J, I ), LDA, WORK(LBWORK*NB+NT*NT+1), 00345 $ K-I+1) 00346 30 CONTINUE 00347 00348 CALL SGEQR2( M-I+1, K-I+1, A( I, I ), LDA, TAU( I ), 00349 $ WORK(LBWORK*NB+NT*NT+1),IINFO ) 00350 00351 ELSE 00352 * 00353 * Use unblocked code to factor the last or only block. 00354 * 00355 CALL SGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), 00356 $ WORK,IINFO ) 00357 00358 END IF 00359 END IF 00360 00361 00362 * 00363 * Apply update to the column M+1:N when N > M 00364 * 00365 IF ( M.LT.N .AND. I.NE.1) THEN 00366 * 00367 * Form the last triangular factor of the block reflector 00368 * H = H(i) H(i+1) . . . H(i+ib-1) 00369 * 00370 IF ( NT .LE. NB ) THEN 00371 CALL SLARFT( 'Forward', 'Columnwise', M-I+1, K-I+1, 00372 $ A( I, I ), LDA, TAU( I ), WORK(I), LBWORK ) 00373 ELSE 00374 CALL SLARFT( 'Forward', 'Columnwise', M-I+1, K-I+1, 00375 $ A( I, I ), LDA, TAU( I ), 00376 $ WORK(LBWORK*NB+1), NT ) 00377 END IF 00378 00379 * 00380 * Apply H' to A(1:M,M+1:N) from the left 00381 * 00382 DO 40 J = 1, K-NX, NB 00383 00384 IB = MIN( K-J+1, NB ) 00385 00386 CALL SLARFB( 'Left', 'Transpose', 'Forward', 00387 $ 'Columnwise', M-J+1, N-M, IB, 00388 $ A( J, J ), LDA, WORK(J), LBWORK, 00389 $ A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1), 00390 $ N-M) 00391 00392 40 CONTINUE 00393 00394 IF ( NT.LE.NB ) THEN 00395 CALL SLARFB( 'Left', 'Transpose', 'Forward', 00396 $ 'Columnwise', M-J+1, N-M, K-J+1, 00397 $ A( J, J ), LDA, WORK(J), LBWORK, 00398 $ A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1), 00399 $ N-M) 00400 ELSE 00401 CALL SLARFB( 'Left', 'Transpose', 'Forward', 00402 $ 'Columnwise', M-J+1, N-M, K-J+1, 00403 $ A( J, J ), LDA, 00404 $ WORK(LBWORK*NB+1), 00405 $ NT, A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1), 00406 $ N-M) 00407 END IF 00408 00409 END IF 00410 00411 WORK( 1 ) = IWS 00412 RETURN 00413 * 00414 * End of SGEQRF 00415 * 00416 END