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

Go to the source code of this file.

Functions/Subroutines

subroutine dgeqp3 (M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
 DGEQP3 More...
 

Function/Subroutine Documentation

subroutine dgeqp3 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEQP3

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

Purpose:
 DGEQP3 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(M,N))
          The scalar factors of the elementary reflectors.
[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. 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 dgeqp3.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  DOUBLE PRECISION 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 dgeqrf, dlaqp2, dlaqps, dormqr, dswap, xerbla
180 * ..
181 * .. External Functions ..
182  INTEGER ilaenv
183  DOUBLE PRECISION dnrm2
184  EXTERNAL ilaenv, dnrm2
185 * ..
186 * .. Intrinsic Functions ..
187  INTRINSIC int, max, min
188 * ..
189 * .. Executable Statements ..
190 *
191 * Test input arguments
192 * ====================
193 *
194  info = 0
195  lquery = ( lwork.EQ.-1 )
196  IF( m.LT.0 ) THEN
197  info = -1
198  ELSE IF( n.LT.0 ) THEN
199  info = -2
200  ELSE IF( lda.LT.max( 1, m ) ) THEN
201  info = -4
202  END IF
203 *
204  IF( info.EQ.0 ) THEN
205  minmn = min( m, n )
206  IF( minmn.EQ.0 ) THEN
207  iws = 1
208  lwkopt = 1
209  ELSE
210  iws = 3*n + 1
211  nb = ilaenv( inb, 'DGEQRF', ' ', m, n, -1, -1 )
212  lwkopt = 2*n + ( n + 1 )*nb
213  END IF
214  work( 1 ) = lwkopt
215 *
216  IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
217  info = -8
218  END IF
219  END IF
220 *
221  IF( info.NE.0 ) THEN
222  CALL xerbla( 'DGEQP3', -info )
223  RETURN
224  ELSE IF( lquery ) THEN
225  RETURN
226  END IF
227 *
228 * Quick return if possible.
229 *
230  IF( minmn.EQ.0 ) THEN
231  RETURN
232  END IF
233 *
234 * Move initial columns up front.
235 *
236  nfxd = 1
237  DO 10 j = 1, n
238  IF( jpvt( j ).NE.0 ) THEN
239  IF( j.NE.nfxd ) THEN
240  CALL dswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
241  jpvt( j ) = jpvt( nfxd )
242  jpvt( nfxd ) = j
243  ELSE
244  jpvt( j ) = j
245  END IF
246  nfxd = nfxd + 1
247  ELSE
248  jpvt( j ) = j
249  END IF
250  10 CONTINUE
251  nfxd = nfxd - 1
252 *
253 * Factorize fixed columns
254 * =======================
255 *
256 * Compute the QR factorization of fixed columns and update
257 * remaining columns.
258 *
259  IF( nfxd.GT.0 ) THEN
260  na = min( m, nfxd )
261 *CC CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
262  CALL dgeqrf( m, na, a, lda, tau, work, lwork, info )
263  iws = max( iws, int( work( 1 ) ) )
264  IF( na.LT.n ) THEN
265 *CC CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
266 *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO )
267  CALL dormqr( 'Left', 'Transpose', m, n-na, na, a, lda, tau,
268  $ a( 1, na+1 ), lda, work, lwork, info )
269  iws = max( iws, int( work( 1 ) ) )
270  END IF
271  END IF
272 *
273 * Factorize free columns
274 * ======================
275 *
276  IF( nfxd.LT.minmn ) THEN
277 *
278  sm = m - nfxd
279  sn = n - nfxd
280  sminmn = minmn - nfxd
281 *
282 * Determine the block size.
283 *
284  nb = ilaenv( inb, 'DGEQRF', ' ', sm, sn, -1, -1 )
285  nbmin = 2
286  nx = 0
287 *
288  IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
289 *
290 * Determine when to cross over from blocked to unblocked code.
291 *
292  nx = max( 0, ilaenv( ixover, 'DGEQRF', ' ', sm, sn, -1,
293  $ -1 ) )
294 *
295 *
296  IF( nx.LT.sminmn ) THEN
297 *
298 * Determine if workspace is large enough for blocked code.
299 *
300  minws = 2*sn + ( sn+1 )*nb
301  iws = max( iws, minws )
302  IF( lwork.LT.minws ) THEN
303 *
304 * Not enough workspace to use optimal NB: Reduce NB and
305 * determine the minimum value of NB.
306 *
307  nb = ( lwork-2*sn ) / ( sn+1 )
308  nbmin = max( 2, ilaenv( inbmin, 'DGEQRF', ' ', sm, sn,
309  $ -1, -1 ) )
310 *
311 *
312  END IF
313  END IF
314  END IF
315 *
316 * Initialize partial column norms. The first N elements of work
317 * store the exact column norms.
318 *
319  DO 20 j = nfxd + 1, n
320  work( j ) = dnrm2( sm, a( nfxd+1, j ), 1 )
321  work( n+j ) = work( j )
322  20 CONTINUE
323 *
324  IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
325  $ ( nx.LT.sminmn ) ) THEN
326 *
327 * Use blocked code initially.
328 *
329  j = nfxd + 1
330 *
331 * Compute factorization: while loop.
332 *
333 *
334  topbmn = minmn - nx
335  30 CONTINUE
336  IF( j.LE.topbmn ) THEN
337  jb = min( nb, topbmn-j+1 )
338 *
339 * Factorize JB columns among columns J:N.
340 *
341  CALL dlaqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
342  $ jpvt( j ), tau( j ), work( j ), work( n+j ),
343  $ work( 2*n+1 ), work( 2*n+jb+1 ), n-j+1 )
344 *
345  j = j + fjb
346  GO TO 30
347  END IF
348  ELSE
349  j = nfxd + 1
350  END IF
351 *
352 * Use unblocked code to factor the last or only block.
353 *
354 *
355  IF( j.LE.minmn )
356  $ CALL dlaqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
357  $ tau( j ), work( j ), work( n+j ),
358  $ work( 2*n+1 ) )
359 *
360  END IF
361 *
362  work( 1 ) = iws
363  RETURN
364 *
365 * End of DGEQP3
366 *
subroutine dlaqp2(M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, WORK)
DLAQP2 computes a QR factorization with column pivoting of the matrix block.
Definition: dlaqp2.f:151
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:138
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:171
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine dlaqps(M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, VN2, AUXV, F, LDF)
DLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BL...
Definition: dlaqps.f:179

Here is the call graph for this function:

Here is the caller graph for this function: