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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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: