LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dgeqrf()

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
December 2016

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 150 of file dgeqrf.f.

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