LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sgeqrf.f
Go to the documentation of this file.
1 C> \brief \b SGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE SGEQRF ( M, N, A, LDA, TAU, WORK, LWORK, INFO )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER INFO, LDA, LWORK, M, N
15 * ..
16 * .. Array Arguments ..
17 * REAL A( LDA, * ), TAU( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 C>\details \b Purpose:
24 C>\verbatim
25 C>
26 C> SGEQRF computes a QR factorization of a real M-by-N matrix A:
27 C> A = Q * R.
28 C>
29 C> This is the left-looking Level 3 BLAS version of the algorithm.
30 C>
31 C>\endverbatim
32 *
33 * Arguments:
34 * ==========
35 *
36 C> \param[in] M
37 C> \verbatim
38 C> M is INTEGER
39 C> The number of rows of the matrix A. M >= 0.
40 C> \endverbatim
41 C>
42 C> \param[in] N
43 C> \verbatim
44 C> N is INTEGER
45 C> The number of columns of the matrix A. N >= 0.
46 C> \endverbatim
47 C>
48 C> \param[in,out] A
49 C> \verbatim
50 C> A is REAL array, dimension (LDA,N)
51 C> On entry, the M-by-N matrix A.
52 C> On exit, the elements on and above the diagonal of the array
53 C> contain the min(M,N)-by-N upper trapezoidal matrix R (R is
54 C> upper triangular if m >= n); the elements below the diagonal,
55 C> with the array TAU, represent the orthogonal matrix Q as a
56 C> product of min(m,n) elementary reflectors (see Further
57 C> Details).
58 C> \endverbatim
59 C>
60 C> \param[in] LDA
61 C> \verbatim
62 C> LDA is INTEGER
63 C> The leading dimension of the array A. LDA >= max(1,M).
64 C> \endverbatim
65 C>
66 C> \param[out] TAU
67 C> \verbatim
68 C> TAU is REAL array, dimension (min(M,N))
69 C> The scalar factors of the elementary reflectors (see Further
70 C> Details).
71 C> \endverbatim
72 C>
73 C> \param[out] WORK
74 C> \verbatim
75 C> WORK is REAL array, dimension (MAX(1,LWORK))
76 C> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
77 C> \endverbatim
78 C>
79 C> \param[in] LWORK
80 C> \verbatim
81 C> LWORK is INTEGER
82 C> \endverbatim
83 C> \verbatim
84 C> The dimension of the array WORK. The dimension can be divided into three parts.
85 C> \endverbatim
86 C> \verbatim
87 C> 1) The part for the triangular factor T. If the very last T is not bigger
88 C> than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
89 C> NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
90 C> \endverbatim
91 C> \verbatim
92 C> 2) The part for the very last T when T is bigger than any of the rest T.
93 C> The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
94 C> where K = min(M,N), NX is calculated by
95 C> NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) )
96 C> \endverbatim
97 C> \verbatim
98 C> 3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
99 C> \endverbatim
100 C> \verbatim
101 C> So LWORK = part1 + part2 + part3
102 C> \endverbatim
103 C> \verbatim
104 C> If LWORK = -1, then a workspace query is assumed; the routine
105 C> only calculates the optimal size of the WORK array, returns
106 C> this value as the first entry of the WORK array, and no error
107 C> message related to LWORK is issued by XERBLA.
108 C> \endverbatim
109 C>
110 C> \param[out] INFO
111 C> \verbatim
112 C> INFO is INTEGER
113 C> = 0: successful exit
114 C> < 0: if INFO = -i, the i-th argument had an illegal value
115 C> \endverbatim
116 C>
117 *
118 * Authors:
119 * ========
120 *
121 C> \author Univ. of Tennessee
122 C> \author Univ. of California Berkeley
123 C> \author Univ. of Colorado Denver
124 C> \author NAG Ltd.
125 *
126 C> \date November 2011
127 *
128 C> \ingroup variantsGEcomputational
129 *
130 * Further Details
131 * ===============
132 C>\details \b Further \b Details
133 C> \verbatim
134 C>
135 C> The matrix Q is represented as a product of elementary reflectors
136 C>
137 C> Q = H(1) H(2) . . . H(k), where k = min(m,n).
138 C>
139 C> Each H(i) has the form
140 C>
141 C> H(i) = I - tau * v * v'
142 C>
143 C> where tau is a real scalar, and v is a real vector with
144 C> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
145 C> and tau in TAU(i).
146 C>
147 C> \endverbatim
148 C>
149 * =====================================================================
150  SUBROUTINE sgeqrf ( M, N, A, LDA, TAU, WORK, LWORK, INFO )
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  REAL 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 sgeqr2, slarfb, slarft, 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, 'SGEQRF', ' ', 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, 'SGEQRF', ' ', 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( 'SGEQRF', -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, 'SGEQRF', ' ', 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 slarfb( '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 sgeqr2( 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 slarft( '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 slarfb( '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 sgeqr2( 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 sgeqr2( 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 slarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
372  $ a( i, i ), lda, tau( i ), work(i), lbwork )
373  ELSE
374  CALL slarft( '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 slarfb( '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 slarfb( '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 slarfb( '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 SGEQRF
415 *
416  END