LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine dgeqrf ( integer M, integer N, double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( * ) TAU, double precision, dimension( * ) WORK, integer LWORK, integer INFO )

DGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.

Purpose:

``` DGEQRF computes a QR factorization of a real M-by-N matrix A:
A = Q * R.

This is the left-looking Level 3 BLAS version of the algorithm.```
Parameters
 [in] M ``` M is INTEGER The number of rows of the matrix A. M >= 0.``` [in] N ``` N is INTEGER The number of columns of the matrix A. N >= 0.``` [in,out] A ``` A is DOUBLE PRECISION array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, the elements on and above the diagonal of the array contain the min(M,N)-by-N upper trapezoidal matrix R (R is upper triangular if m >= n); the elements below the diagonal, with the array TAU, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors (see Further Details).``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).``` [out] TAU ``` TAU is DOUBLE PRECISION array, dimension (min(M,N)) The scalar factors of the elementary reflectors (see Further Details).``` [out] WORK ``` WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.``` [in] LWORK ` LWORK is INTEGER` ` The dimension of the array WORK. The dimension can be divided into three parts.` ``` 1) The part for the triangular factor T. If the very last T is not bigger than any of the rest, then this part is NB x ceiling(K/NB), otherwise, NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T ``` ``` 2) The part for the very last T when T is bigger than any of the rest T. The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB, where K = min(M,N), NX is calculated by NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )``` ` 3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)` ` So LWORK = part1 + part2 + part3` ``` If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Date
November 2011

Further Details

```  The matrix Q is represented as a product of elementary reflectors

Q = H(1) H(2) . . . H(k), where k = min(m,n).

Each H(i) has the form

H(i) = I - tau * v * v'

where tau is a real scalar, and v is a real vector with
v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
and tau in TAU(i).```

Definition at line 151 of file dgeqrf.f.

151 *
152 * -- LAPACK computational routine (version 3.1) --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155 * November 2011
156 *
157 * .. Scalar Arguments ..
158  INTEGER info, lda, lwork, m, n
159 * ..
160 * .. Array Arguments ..
161  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Local Scalars ..
167  LOGICAL lquery
168  INTEGER i, ib, iinfo, iws, j, k, lwkopt, nb,
169  \$ nbmin, nx, lbwork, nt, llwork
170 * ..
171 * .. External Subroutines ..
172  EXTERNAL dgeqr2, dlarfb, dlarft, xerbla
173 * ..
174 * .. Intrinsic Functions ..
175  INTRINSIC max, min
176 * ..
177 * .. External Functions ..
178  INTEGER ilaenv
179  REAL sceil
180  EXTERNAL ilaenv, sceil
181 * ..
182 * .. Executable Statements ..
183
184  info = 0
185  nbmin = 2
186  nx = 0
187  iws = n
188  k = min( m, n )
189  nb = ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
190
191  IF( nb.GT.1 .AND. nb.LT.k ) THEN
192 *
193 * Determine when to cross over from blocked to unblocked code.
194 *
195  nx = max( 0, ilaenv( 3, 'DGEQRF', ' ', m, n, -1, -1 ) )
196  END IF
197 *
198 * Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
199 *
200 * NB=3 2NB=6 K=10
201 * | | |
202 * 1--2--3--4--5--6--7--8--9--10
203 * | \________/
204 * K-NX=5 NT=4
205 *
206 * So here 4 x 4 is the last T stored in the workspace
207 *
208  nt = k-sceil(REAL(k-nx)/REAL(nb))*nb
209
210 *
211 * optimal workspace = space for dlarfb + space for normal T's + space for the last T
212 *
213  llwork = max(max((n-m)*k, (n-m)*nb), max(k*nb, nb*nb))
214  llwork = sceil(REAL(llwork)/REAL(nb))
215
216  IF ( nt.GT.nb ) THEN
217
218  lbwork = k-nt
219 *
220 * Optimal workspace for dlarfb = MAX(1,N)*NT
221 *
222  lwkopt = (lbwork+llwork)*nb
223  work( 1 ) = (lwkopt+nt*nt)
224
225  ELSE
226
227  lbwork = sceil(REAL(k)/REAL(nb))*nb
228  lwkopt = (lbwork+llwork-nb)*nb
229  work( 1 ) = lwkopt
230
231  END IF
232
233 *
234 * Test the input arguments
235 *
236  lquery = ( lwork.EQ.-1 )
237  IF( m.LT.0 ) THEN
238  info = -1
239  ELSE IF( n.LT.0 ) THEN
240  info = -2
241  ELSE IF( lda.LT.max( 1, m ) ) THEN
242  info = -4
243  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
244  info = -7
245  END IF
246  IF( info.NE.0 ) THEN
247  CALL xerbla( 'DGEQRF', -info )
248  RETURN
249  ELSE IF( lquery ) THEN
250  RETURN
251  END IF
252 *
253 * Quick return if possible
254 *
255  IF( k.EQ.0 ) THEN
256  work( 1 ) = 1
257  RETURN
258  END IF
259 *
260  IF( nb.GT.1 .AND. nb.LT.k ) THEN
261
262  IF( nx.LT.k ) THEN
263 *
264 * Determine if workspace is large enough for blocked code.
265 *
266  IF ( nt.LE.nb ) THEN
267  iws = (lbwork+llwork-nb)*nb
268  ELSE
269  iws = (lbwork+llwork)*nb+nt*nt
270  END IF
271
272  IF( lwork.LT.iws ) THEN
273 *
274 * Not enough workspace to use optimal NB: reduce NB and
275 * determine the minimum value of NB.
276 *
277  IF ( nt.LE.nb ) THEN
278  nb = lwork / (llwork+(lbwork-nb))
279  ELSE
280  nb = (lwork-nt*nt)/(lbwork+llwork)
281  END IF
282
283  nbmin = max( 2, ilaenv( 2, 'DGEQRF', ' ', m, n, -1,
284  \$ -1 ) )
285  END IF
286  END IF
287  END IF
288 *
289  IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
290 *
291 * Use blocked code initially
292 *
293  DO 10 i = 1, k - nx, nb
294  ib = min( k-i+1, nb )
295 *
296 * Update the current column using old T's
297 *
298  DO 20 j = 1, i - nb, nb
299 *
300 * Apply H' to A(J:M,I:I+IB-1) from the left
301 *
302  CALL dlarfb( 'Left', 'Transpose', 'Forward',
303  \$ 'Columnwise', m-j+1, ib, nb,
304  \$ a( j, j ), lda, work(j), lbwork,
305  \$ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
306  \$ ib)
307
308 20 CONTINUE
309 *
310 * Compute the QR factorization of the current block
311 * A(I:M,I:I+IB-1)
312 *
313  CALL dgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
314  \$ work(lbwork*nb+nt*nt+1), iinfo )
315
316  IF( i+ib.LE.n ) THEN
317 *
318 * Form the triangular factor of the block reflector
319 * H = H(i) H(i+1) . . . H(i+ib-1)
320 *
321  CALL dlarft( 'Forward', 'Columnwise', m-i+1, ib,
322  \$ a( i, i ), lda, tau( i ),
323  \$ work(i), lbwork )
324 *
325  END IF
326  10 CONTINUE
327  ELSE
328  i = 1
329  END IF
330 *
331 * Use unblocked code to factor the last or only block.
332 *
333  IF( i.LE.k ) THEN
334
335  IF ( i .NE. 1 ) THEN
336
337  DO 30 j = 1, i - nb, nb
338 *
339 * Apply H' to A(J:M,I:K) from the left
340 *
341  CALL dlarfb( 'Left', 'Transpose', 'Forward',
342  \$ 'Columnwise', m-j+1, k-i+1, nb,
343  \$ a( j, j ), lda, work(j), lbwork,
344  \$ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
345  \$ k-i+1)
346 30 CONTINUE
347
348  CALL dgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
349  \$ work(lbwork*nb+nt*nt+1),iinfo )
350
351  ELSE
352 *
353 * Use unblocked code to factor the last or only block.
354 *
355  CALL dgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
356  \$ work,iinfo )
357
358  END IF
359  END IF
360
361
362 *
363 * Apply update to the column M+1:N when N > M
364 *
365  IF ( m.LT.n .AND. i.NE.1) THEN
366 *
367 * Form the last triangular factor of the block reflector
368 * H = H(i) H(i+1) . . . H(i+ib-1)
369 *
370  IF ( nt .LE. nb ) THEN
371  CALL dlarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
372  \$ a( i, i ), lda, tau( i ), work(i), lbwork )
373  ELSE
374  CALL dlarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
375  \$ a( i, i ), lda, tau( i ),
376  \$ work(lbwork*nb+1), nt )
377  END IF
378
379 *
380 * Apply H' to A(1:M,M+1:N) from the left
381 *
382  DO 40 j = 1, k-nx, nb
383
384  ib = min( k-j+1, nb )
385
386  CALL dlarfb( 'Left', 'Transpose', 'Forward',
387  \$ 'Columnwise', m-j+1, n-m, ib,
388  \$ a( j, j ), lda, work(j), lbwork,
389  \$ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
390  \$ n-m)
391
392 40 CONTINUE
393
394  IF ( nt.LE.nb ) THEN
395  CALL dlarfb( 'Left', 'Transpose', 'Forward',
396  \$ 'Columnwise', m-j+1, n-m, k-j+1,
397  \$ a( j, j ), lda, work(j), lbwork,
398  \$ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
399  \$ n-m)
400  ELSE
401  CALL dlarfb( 'Left', 'Transpose', 'Forward',
402  \$ 'Columnwise', m-j+1, n-m, k-j+1,
403  \$ a( j, j ), lda,
404  \$ work(lbwork*nb+1),
405  \$ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
406  \$ n-m)
407  END IF
408
409  END IF
410
411  work( 1 ) = iws
412  RETURN
413 *
414 * End of DGEQRF
415 *
subroutine dlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition: dlarfb.f:197
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgeqr2.f:123
real function sceil(A)
SCEIL
Definition: sceil.f:60
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: dlarft.f:165
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function: