LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
sgeqp3.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine sgeqp3 (M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
 SGEQP3 More...
 

Function/Subroutine Documentation

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.
Date
September 2012
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 153 of file sgeqp3.f.

153 *
154 * -- LAPACK computational routine (version 3.4.2) --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 * September 2012
158 *
159 * .. Scalar Arguments ..
160  INTEGER info, lda, lwork, m, n
161 * ..
162 * .. Array Arguments ..
163  INTEGER jpvt( * )
164  REAL a( lda, * ), tau( * ), work( * )
165 * ..
166 *
167 * =====================================================================
168 *
169 * .. Parameters ..
170  INTEGER inb, inbmin, ixover
171  parameter( inb = 1, inbmin = 2, ixover = 3 )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL lquery
175  INTEGER fjb, iws, j, jb, lwkopt, minmn, minws, na, nb,
176  $ nbmin, nfxd, nx, sm, sminmn, sn, topbmn
177 * ..
178 * .. External Subroutines ..
179  EXTERNAL sgeqrf, slaqp2, slaqps, sormqr, sswap, xerbla
180 * ..
181 * .. External Functions ..
182  INTEGER ilaenv
183  REAL snrm2
184  EXTERNAL ilaenv, snrm2
185 * ..
186 * .. Intrinsic Functions ..
187  INTRINSIC int, max, min
188 * ..
189 * .. Executable Statements ..
190 *
191  info = 0
192  lquery = ( lwork.EQ.-1 )
193  IF( m.LT.0 ) THEN
194  info = -1
195  ELSE IF( n.LT.0 ) THEN
196  info = -2
197  ELSE IF( lda.LT.max( 1, m ) ) THEN
198  info = -4
199  END IF
200 *
201  IF( info.EQ.0 ) THEN
202  minmn = min( m, n )
203  IF( minmn.EQ.0 ) THEN
204  iws = 1
205  lwkopt = 1
206  ELSE
207  iws = 3*n + 1
208  nb = ilaenv( inb, 'SGEQRF', ' ', m, n, -1, -1 )
209  lwkopt = 2*n + ( n + 1 )*nb
210  END IF
211  work( 1 ) = lwkopt
212 *
213  IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
214  info = -8
215  END IF
216  END IF
217 *
218  IF( info.NE.0 ) THEN
219  CALL xerbla( 'SGEQP3', -info )
220  RETURN
221  ELSE IF( lquery ) THEN
222  RETURN
223  END IF
224 *
225 * Quick return if possible.
226 *
227  IF( minmn.EQ.0 ) THEN
228  RETURN
229  END IF
230 *
231 * Move initial columns up front.
232 *
233  nfxd = 1
234  DO 10 j = 1, n
235  IF( jpvt( j ).NE.0 ) THEN
236  IF( j.NE.nfxd ) THEN
237  CALL sswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
238  jpvt( j ) = jpvt( nfxd )
239  jpvt( nfxd ) = j
240  ELSE
241  jpvt( j ) = j
242  END IF
243  nfxd = nfxd + 1
244  ELSE
245  jpvt( j ) = j
246  END IF
247  10 CONTINUE
248  nfxd = nfxd - 1
249 *
250 * Factorize fixed columns
251 * =======================
252 *
253 * Compute the QR factorization of fixed columns and update
254 * remaining columns.
255 *
256  IF( nfxd.GT.0 ) THEN
257  na = min( m, nfxd )
258 *CC CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
259  CALL sgeqrf( m, na, a, lda, tau, work, lwork, info )
260  iws = max( iws, int( work( 1 ) ) )
261  IF( na.LT.n ) THEN
262 *CC CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
263 *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO )
264  CALL sormqr( 'Left', 'Transpose', m, n-na, na, a, lda, tau,
265  $ a( 1, na+1 ), lda, work, lwork, info )
266  iws = max( iws, int( work( 1 ) ) )
267  END IF
268  END IF
269 *
270 * Factorize free columns
271 * ======================
272 *
273  IF( nfxd.LT.minmn ) THEN
274 *
275  sm = m - nfxd
276  sn = n - nfxd
277  sminmn = minmn - nfxd
278 *
279 * Determine the block size.
280 *
281  nb = ilaenv( inb, 'SGEQRF', ' ', sm, sn, -1, -1 )
282  nbmin = 2
283  nx = 0
284 *
285  IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
286 *
287 * Determine when to cross over from blocked to unblocked code.
288 *
289  nx = max( 0, ilaenv( ixover, 'SGEQRF', ' ', sm, sn, -1,
290  $ -1 ) )
291 *
292 *
293  IF( nx.LT.sminmn ) THEN
294 *
295 * Determine if workspace is large enough for blocked code.
296 *
297  minws = 2*sn + ( sn+1 )*nb
298  iws = max( iws, minws )
299  IF( lwork.LT.minws ) THEN
300 *
301 * Not enough workspace to use optimal NB: Reduce NB and
302 * determine the minimum value of NB.
303 *
304  nb = ( lwork-2*sn ) / ( sn+1 )
305  nbmin = max( 2, ilaenv( inbmin, 'SGEQRF', ' ', sm, sn,
306  $ -1, -1 ) )
307 *
308 *
309  END IF
310  END IF
311  END IF
312 *
313 * Initialize partial column norms. The first N elements of work
314 * store the exact column norms.
315 *
316  DO 20 j = nfxd + 1, n
317  work( j ) = snrm2( sm, a( nfxd+1, j ), 1 )
318  work( n+j ) = work( j )
319  20 CONTINUE
320 *
321  IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
322  $ ( nx.LT.sminmn ) ) THEN
323 *
324 * Use blocked code initially.
325 *
326  j = nfxd + 1
327 *
328 * Compute factorization: while loop.
329 *
330 *
331  topbmn = minmn - nx
332  30 CONTINUE
333  IF( j.LE.topbmn ) THEN
334  jb = min( nb, topbmn-j+1 )
335 *
336 * Factorize JB columns among columns J:N.
337 *
338  CALL slaqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
339  $ jpvt( j ), tau( j ), work( j ), work( n+j ),
340  $ work( 2*n+1 ), work( 2*n+jb+1 ), n-j+1 )
341 *
342  j = j + fjb
343  GO TO 30
344  END IF
345  ELSE
346  j = nfxd + 1
347  END IF
348 *
349 * Use unblocked code to factor the last or only block.
350 *
351 *
352  IF( j.LE.minmn )
353  $ CALL slaqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
354  $ tau( j ), work( j ), work( n+j ),
355  $ work( 2*n+1 ) )
356 *
357  END IF
358 *
359  work( 1 ) = iws
360  RETURN
361 *
362 * End of SGEQP3
363 *
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:172
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:151
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:180
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:138
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83

Here is the call graph for this function:

Here is the caller graph for this function: