LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cgeqrf()

subroutine cgeqrf ( integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( * )  TAU,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

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

Purpose:

 CGEQRF 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 COMPLEX 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 COMPLEX array, dimension (min(M,N))
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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. LWORK >= 1 if MIN(M,N) = 0,
          otherwise 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, 'CGEQRF', ' ', 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 151 of file cgeqrf.f.

152*
153* -- LAPACK computational routine --
154* -- LAPACK is a software package provided by Univ. of Tennessee, --
155* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156*
157* .. Scalar Arguments ..
158 INTEGER INFO, LDA, LWORK, M, N
159* ..
160* .. Array Arguments ..
161 COMPLEX 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 cgeqr2, clarfb, clarft, 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, 'CGEQRF', ' ', 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, 'CGEQRF', ' ', 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( k.EQ.0 ) THEN
217
218 lbwork = 0
219 lwkopt = 1
220 work( 1 ) = lwkopt
221
222 ELSE IF ( nt.GT.nb ) THEN
223
224 lbwork = k-nt
225*
226* Optimal workspace for dlarfb = MAX(1,N)*NT
227*
228 lwkopt = (lbwork+llwork)*nb
229 work( 1 ) = (lwkopt+nt*nt)
230
231 ELSE
232
233 lbwork = sceil(real(k)/real(nb))*nb
234 lwkopt = (lbwork+llwork-nb)*nb
235 work( 1 ) = lwkopt
236
237 END IF
238
239*
240* Test the input arguments
241*
242 lquery = ( lwork.EQ.-1 )
243 IF( m.LT.0 ) THEN
244 info = -1
245 ELSE IF( n.LT.0 ) THEN
246 info = -2
247 ELSE IF( lda.LT.max( 1, m ) ) THEN
248 info = -4
249 ELSE IF ( .NOT.lquery ) THEN
250 IF( lwork.LE.0 .OR. ( m.GT.0 .AND. lwork.LT.max( 1, n ) ) )
251 $ info = -7
252 END IF
253 IF( info.NE.0 ) THEN
254 CALL xerbla( 'CGEQRF', -info )
255 RETURN
256 ELSE IF( lquery ) THEN
257 RETURN
258 END IF
259*
260* Quick return if possible
261*
262 IF( k.EQ.0 ) THEN
263 RETURN
264 END IF
265*
266 IF( nb.GT.1 .AND. nb.LT.k ) THEN
267
268 IF( nx.LT.k ) THEN
269*
270* Determine if workspace is large enough for blocked code.
271*
272 IF ( nt.LE.nb ) THEN
273 iws = (lbwork+llwork-nb)*nb
274 ELSE
275 iws = (lbwork+llwork)*nb+nt*nt
276 END IF
277
278 IF( lwork.LT.iws ) THEN
279*
280* Not enough workspace to use optimal NB: reduce NB and
281* determine the minimum value of NB.
282*
283 IF ( nt.LE.nb ) THEN
284 nb = lwork / (llwork+(lbwork-nb))
285 ELSE
286 nb = (lwork-nt*nt)/(lbwork+llwork)
287 END IF
288
289 nbmin = max( 2, ilaenv( 2, 'CGEQRF', ' ', m, n, -1,
290 $ -1 ) )
291 END IF
292 END IF
293 END IF
294*
295 IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
296*
297* Use blocked code initially
298*
299 DO 10 i = 1, k - nx, nb
300 ib = min( k-i+1, nb )
301*
302* Update the current column using old T's
303*
304 DO 20 j = 1, i - nb, nb
305*
306* Apply H' to A(J:M,I:I+IB-1) from the left
307*
308 CALL clarfb( 'Left', 'Transpose', 'Forward',
309 $ 'Columnwise', m-j+1, ib, nb,
310 $ a( j, j ), lda, work(j), lbwork,
311 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
312 $ ib)
313
31420 CONTINUE
315*
316* Compute the QR factorization of the current block
317* A(I:M,I:I+IB-1)
318*
319 CALL cgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
320 $ work(lbwork*nb+nt*nt+1), iinfo )
321
322 IF( i+ib.LE.n ) THEN
323*
324* Form the triangular factor of the block reflector
325* H = H(i) H(i+1) . . . H(i+ib-1)
326*
327 CALL clarft( 'Forward', 'Columnwise', m-i+1, ib,
328 $ a( i, i ), lda, tau( i ),
329 $ work(i), lbwork )
330*
331 END IF
332 10 CONTINUE
333 ELSE
334 i = 1
335 END IF
336*
337* Use unblocked code to factor the last or only block.
338*
339 IF( i.LE.k ) THEN
340
341 IF ( i .NE. 1 ) THEN
342
343 DO 30 j = 1, i - nb, nb
344*
345* Apply H' to A(J:M,I:K) from the left
346*
347 CALL clarfb( 'Left', 'Transpose', 'Forward',
348 $ 'Columnwise', m-j+1, k-i+1, nb,
349 $ a( j, j ), lda, work(j), lbwork,
350 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
351 $ k-i+1)
35230 CONTINUE
353
354 CALL cgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
355 $ work(lbwork*nb+nt*nt+1),iinfo )
356
357 ELSE
358*
359* Use unblocked code to factor the last or only block.
360*
361 CALL cgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
362 $ work,iinfo )
363
364 END IF
365 END IF
366
367
368*
369* Apply update to the column M+1:N when N > M
370*
371 IF ( m.LT.n .AND. i.NE.1) THEN
372*
373* Form the last triangular factor of the block reflector
374* H = H(i) H(i+1) . . . H(i+ib-1)
375*
376 IF ( nt .LE. nb ) THEN
377 CALL clarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
378 $ a( i, i ), lda, tau( i ), work(i), lbwork )
379 ELSE
380 CALL clarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
381 $ a( i, i ), lda, tau( i ),
382 $ work(lbwork*nb+1), nt )
383 END IF
384
385*
386* Apply H' to A(1:M,M+1:N) from the left
387*
388 DO 40 j = 1, k-nx, nb
389
390 ib = min( k-j+1, nb )
391
392 CALL clarfb( 'Left', 'Transpose', 'Forward',
393 $ 'Columnwise', m-j+1, n-m, ib,
394 $ a( j, j ), lda, work(j), lbwork,
395 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
396 $ n-m)
397
39840 CONTINUE
399
400 IF ( nt.LE.nb ) THEN
401 CALL clarfb( 'Left', 'Transpose', 'Forward',
402 $ 'Columnwise', m-j+1, n-m, k-j+1,
403 $ a( j, j ), lda, work(j), lbwork,
404 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
405 $ n-m)
406 ELSE
407 CALL clarfb( 'Left', 'Transpose', 'Forward',
408 $ 'Columnwise', m-j+1, n-m, k-j+1,
409 $ a( j, j ), lda,
410 $ work(lbwork*nb+1),
411 $ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
412 $ n-m)
413 END IF
414
415 END IF
416
417 work( 1 ) = iws
418 RETURN
419*
420* End of CGEQRF
421*
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 cgeqr2(M, N, A, LDA, TAU, WORK, INFO)
CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition: cgeqr2.f:130
subroutine clarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
Definition: clarfb.f:197
subroutine clarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: clarft.f:163
real function sceil(A)
SCEIL
Definition: sceil.f:60
Here is the call graph for this function: