LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sgeqp3()

subroutine sgeqp3 ( integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SGEQP3

Download SGEQP3 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SGEQP3 computes a QR factorization with column pivoting of a
 matrix A:  A*P = Q*R  using Level 3 BLAS.
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 REAL array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, the upper triangle of the array contains the
          min(M,N)-by-N upper trapezoidal matrix R; the elements below
          the diagonal, together with the array TAU, represent the
          orthogonal matrix Q as a product of min(M,N) elementary
          reflectors.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in,out]JPVT
          JPVT is INTEGER array, dimension (N)
          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
          to the front of A*P (a leading column); if JPVT(J)=0,
          the J-th column of A is a free column.
          On exit, if JPVT(J)=K, then the J-th column of A*P was the
          the K-th column of A.
[out]TAU
          TAU is REAL array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[out]WORK
          WORK is REAL 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 >= 3*N+1.
          For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB
          is the optimal blocksize.

          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.
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**T

  where tau is a real scalar, and v is a real/complex 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).
Contributors:
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA

Definition at line 150 of file sgeqp3.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  INTEGER JPVT( * )
161  REAL A( LDA, * ), TAU( * ), WORK( * )
162 * ..
163 *
164 * =====================================================================
165 *
166 * .. Parameters ..
167  INTEGER INB, INBMIN, IXOVER
168  parameter( inb = 1, inbmin = 2, ixover = 3 )
169 * ..
170 * .. Local Scalars ..
171  LOGICAL LQUERY
172  INTEGER FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
173  $ NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL sgeqrf, slaqp2, slaqps, sormqr, sswap, xerbla
177 * ..
178 * .. External Functions ..
179  INTEGER ILAENV
180  REAL SNRM2
181  EXTERNAL ilaenv, snrm2
182 * ..
183 * .. Intrinsic Functions ..
184  INTRINSIC int, max, min
185 * Test input arguments
186 * ====================
187 *
188  info = 0
189  lquery = ( lwork.EQ.-1 )
190  IF( m.LT.0 ) THEN
191  info = -1
192  ELSE IF( n.LT.0 ) THEN
193  info = -2
194  ELSE IF( lda.LT.max( 1, m ) ) THEN
195  info = -4
196  END IF
197 *
198  IF( info.EQ.0 ) THEN
199  minmn = min( m, n )
200  IF( minmn.EQ.0 ) THEN
201  iws = 1
202  lwkopt = 1
203  ELSE
204  iws = 3*n + 1
205  nb = ilaenv( inb, 'SGEQRF', ' ', m, n, -1, -1 )
206  lwkopt = 2*n + ( n + 1 )*nb
207  END IF
208  work( 1 ) = lwkopt
209 *
210  IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
211  info = -8
212  END IF
213  END IF
214 *
215  IF( info.NE.0 ) THEN
216  CALL xerbla( 'SGEQP3', -info )
217  RETURN
218  ELSE IF( lquery ) THEN
219  RETURN
220  END IF
221 *
222 * Move initial columns up front.
223 *
224  nfxd = 1
225  DO 10 j = 1, n
226  IF( jpvt( j ).NE.0 ) THEN
227  IF( j.NE.nfxd ) THEN
228  CALL sswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
229  jpvt( j ) = jpvt( nfxd )
230  jpvt( nfxd ) = j
231  ELSE
232  jpvt( j ) = j
233  END IF
234  nfxd = nfxd + 1
235  ELSE
236  jpvt( j ) = j
237  END IF
238  10 CONTINUE
239  nfxd = nfxd - 1
240 *
241 * Factorize fixed columns
242 * =======================
243 *
244 * Compute the QR factorization of fixed columns and update
245 * remaining columns.
246 *
247  IF( nfxd.GT.0 ) THEN
248  na = min( m, nfxd )
249 *CC CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
250  CALL sgeqrf( m, na, a, lda, tau, work, lwork, info )
251  iws = max( iws, int( work( 1 ) ) )
252  IF( na.LT.n ) THEN
253 *CC CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
254 *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO )
255  CALL sormqr( 'Left', 'Transpose', m, n-na, na, a, lda, tau,
256  $ a( 1, na+1 ), lda, work, lwork, info )
257  iws = max( iws, int( work( 1 ) ) )
258  END IF
259  END IF
260 *
261 * Factorize free columns
262 * ======================
263 *
264  IF( nfxd.LT.minmn ) THEN
265 *
266  sm = m - nfxd
267  sn = n - nfxd
268  sminmn = minmn - nfxd
269 *
270 * Determine the block size.
271 *
272  nb = ilaenv( inb, 'SGEQRF', ' ', sm, sn, -1, -1 )
273  nbmin = 2
274  nx = 0
275 *
276  IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
277 *
278 * Determine when to cross over from blocked to unblocked code.
279 *
280  nx = max( 0, ilaenv( ixover, 'SGEQRF', ' ', sm, sn, -1,
281  $ -1 ) )
282 *
283 *
284  IF( nx.LT.sminmn ) THEN
285 *
286 * Determine if workspace is large enough for blocked code.
287 *
288  minws = 2*sn + ( sn+1 )*nb
289  iws = max( iws, minws )
290  IF( lwork.LT.minws ) THEN
291 *
292 * Not enough workspace to use optimal NB: Reduce NB and
293 * determine the minimum value of NB.
294 *
295  nb = ( lwork-2*sn ) / ( sn+1 )
296  nbmin = max( 2, ilaenv( inbmin, 'SGEQRF', ' ', sm, sn,
297  $ -1, -1 ) )
298 *
299 *
300  END IF
301  END IF
302  END IF
303 *
304 * Initialize partial column norms. The first N elements of work
305 * store the exact column norms.
306 *
307  DO 20 j = nfxd + 1, n
308  work( j ) = snrm2( sm, a( nfxd+1, j ), 1 )
309  work( n+j ) = work( j )
310  20 CONTINUE
311 *
312  IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
313  $ ( nx.LT.sminmn ) ) THEN
314 *
315 * Use blocked code initially.
316 *
317  j = nfxd + 1
318 *
319 * Compute factorization: while loop.
320 *
321 *
322  topbmn = minmn - nx
323  30 CONTINUE
324  IF( j.LE.topbmn ) THEN
325  jb = min( nb, topbmn-j+1 )
326 *
327 * Factorize JB columns among columns J:N.
328 *
329  CALL slaqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
330  $ jpvt( j ), tau( j ), work( j ), work( n+j ),
331  $ work( 2*n+1 ), work( 2*n+jb+1 ), n-j+1 )
332 *
333  j = j + fjb
334  GO TO 30
335  END IF
336  ELSE
337  j = nfxd + 1
338  END IF
339 *
340 * Use unblocked code to factor the last or only block.
341 *
342 *
343  IF( j.LE.minmn )
344  $ CALL slaqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
345  $ tau( j ), work( j ), work( n+j ),
346  $ work( 2*n+1 ) )
347 *
348  END IF
349 *
350  work( 1 ) = iws
351  RETURN
352 *
353 * End of SGEQP3
354 *
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:145
subroutine slaqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
SLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition: slaqp2.f:149
subroutine slaqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
SLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BL...
Definition: slaqps.f:178
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: