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