150 SUBROUTINE dgeqrf ( M, N, A, LDA, TAU, WORK, LWORK, INFO )
158 INTEGER INFO, LDA, LWORK, M, N
161 DOUBLE PRECISION A( lda, * ), TAU( * ), WORK( * )
168 INTEGER I, IB, IINFO, IWS, J, K, LWKOPT, NB,
169 $ nbmin, nx, lbwork, nt, llwork
180 EXTERNAL ilaenv, sceil
189 nb = ilaenv( 1,
'DGEQRF',
' ', m, n, -1, -1 )
191 IF( nb.GT.1 .AND. nb.LT.k )
THEN
195 nx = max( 0, ilaenv( 3,
'DGEQRF',
' ', m, n, -1, -1 ) )
208 nt = k-sceil(
REAL(k-nx)/
REAL(nb))*nb
213 llwork = max(max((n-m)*k, (n-m)*nb), max(k*nb, nb*nb))
214 llwork = sceil(
REAL(llwork)/
REAL(nb))
222 lwkopt = (lbwork+llwork)*nb
223 work( 1 ) = (lwkopt+nt*nt)
227 lbwork = sceil(
REAL(k)/
REAL(nb))*nb
228 lwkopt = (lbwork+llwork-nb)*nb
236 lquery = ( lwork.EQ.-1 )
239 ELSE IF( n.LT.0 )
THEN
241 ELSE IF( lda.LT.max( 1, m ) )
THEN
243 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery )
THEN
247 CALL xerbla(
'DGEQRF', -info )
249 ELSE IF( lquery )
THEN
260 IF( nb.GT.1 .AND. nb.LT.k )
THEN
267 iws = (lbwork+llwork-nb)*nb
269 iws = (lbwork+llwork)*nb+nt*nt
272 IF( lwork.LT.iws )
THEN
278 nb = lwork / (llwork+(lbwork-nb))
280 nb = (lwork-nt*nt)/(lbwork+llwork)
283 nbmin = max( 2, ilaenv( 2,
'DGEQRF',
' ', m, n, -1,
289 IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k )
THEN
293 DO 10 i = 1, k - nx, nb
294 ib = min( k-i+1, nb )
298 DO 20 j = 1, i - nb, nb
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),
313 CALL dgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
314 $ work(lbwork*nb+nt*nt+1), iinfo )
321 CALL dlarft(
'Forward',
'Columnwise', m-i+1, ib,
322 $ a( i, i ), lda, tau( i ),
337 DO 30 j = 1, i - nb, nb
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),
348 CALL dgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
349 $ work(lbwork*nb+nt*nt+1),iinfo )
355 CALL dgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
365 IF ( m.LT.n .AND. i.NE.1)
THEN
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 )
374 CALL dlarft(
'Forward',
'Columnwise', m-i+1, k-i+1,
375 $ a( i, i ), lda, tau( i ),
376 $ work(lbwork*nb+1), nt )
382 DO 40 j = 1, k-nx, nb
384 ib = min( k-j+1, nb )
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),
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),
401 CALL dlarfb(
'Left',
'Transpose',
'Forward',
402 $
'Columnwise', m-j+1, n-m, k-j+1,
405 $ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
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.
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...