LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dgeqp3()

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
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**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.7.0) --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 * December 2016
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 * Move initial columns up front.
229 *
230  nfxd = 1
231  DO 10 j = 1, n
232  IF( jpvt( j ).NE.0 ) THEN
233  IF( j.NE.nfxd ) THEN
234  CALL dswap( m, a( 1, j ), 1, a( 1, nfxd ), 1 )
235  jpvt( j ) = jpvt( nfxd )
236  jpvt( nfxd ) = j
237  ELSE
238  jpvt( j ) = j
239  END IF
240  nfxd = nfxd + 1
241  ELSE
242  jpvt( j ) = j
243  END IF
244  10 CONTINUE
245  nfxd = nfxd - 1
246 *
247 * Factorize fixed columns
248 * =======================
249 *
250 * Compute the QR factorization of fixed columns and update
251 * remaining columns.
252 *
253  IF( nfxd.GT.0 ) THEN
254  na = min( m, nfxd )
255 *CC CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
256  CALL dgeqrf( m, na, a, lda, tau, work, lwork, info )
257  iws = max( iws, int( work( 1 ) ) )
258  IF( na.LT.n ) THEN
259 *CC CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
260 *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO )
261  CALL dormqr( 'Left', 'Transpose', m, n-na, na, a, lda, tau,
262  $ a( 1, na+1 ), lda, work, lwork, info )
263  iws = max( iws, int( work( 1 ) ) )
264  END IF
265  END IF
266 *
267 * Factorize free columns
268 * ======================
269 *
270  IF( nfxd.LT.minmn ) THEN
271 *
272  sm = m - nfxd
273  sn = n - nfxd
274  sminmn = minmn - nfxd
275 *
276 * Determine the block size.
277 *
278  nb = ilaenv( inb, 'DGEQRF', ' ', sm, sn, -1, -1 )
279  nbmin = 2
280  nx = 0
281 *
282  IF( ( nb.GT.1 ) .AND. ( nb.LT.sminmn ) ) THEN
283 *
284 * Determine when to cross over from blocked to unblocked code.
285 *
286  nx = max( 0, ilaenv( ixover, 'DGEQRF', ' ', sm, sn, -1,
287  $ -1 ) )
288 *
289 *
290  IF( nx.LT.sminmn ) THEN
291 *
292 * Determine if workspace is large enough for blocked code.
293 *
294  minws = 2*sn + ( sn+1 )*nb
295  iws = max( iws, minws )
296  IF( lwork.LT.minws ) THEN
297 *
298 * Not enough workspace to use optimal NB: Reduce NB and
299 * determine the minimum value of NB.
300 *
301  nb = ( lwork-2*sn ) / ( sn+1 )
302  nbmin = max( 2, ilaenv( inbmin, 'DGEQRF', ' ', sm, sn,
303  $ -1, -1 ) )
304 *
305 *
306  END IF
307  END IF
308  END IF
309 *
310 * Initialize partial column norms. The first N elements of work
311 * store the exact column norms.
312 *
313  DO 20 j = nfxd + 1, n
314  work( j ) = dnrm2( sm, a( nfxd+1, j ), 1 )
315  work( n+j ) = work( j )
316  20 CONTINUE
317 *
318  IF( ( nb.GE.nbmin ) .AND. ( nb.LT.sminmn ) .AND.
319  $ ( nx.LT.sminmn ) ) THEN
320 *
321 * Use blocked code initially.
322 *
323  j = nfxd + 1
324 *
325 * Compute factorization: while loop.
326 *
327 *
328  topbmn = minmn - nx
329  30 CONTINUE
330  IF( j.LE.topbmn ) THEN
331  jb = min( nb, topbmn-j+1 )
332 *
333 * Factorize JB columns among columns J:N.
334 *
335  CALL dlaqps( m, n-j+1, j-1, jb, fjb, a( 1, j ), lda,
336  $ jpvt( j ), tau( j ), work( j ), work( n+j ),
337  $ work( 2*n+1 ), work( 2*n+jb+1 ), n-j+1 )
338 *
339  j = j + fjb
340  GO TO 30
341  END IF
342  ELSE
343  j = nfxd + 1
344  END IF
345 *
346 * Use unblocked code to factor the last or only block.
347 *
348 *
349  IF( j.LE.minmn )
350  $ CALL dlaqp2( m, n-j+1, j-1, a( 1, j ), lda, jpvt( j ),
351  $ tau( j ), work( j ), work( n+j ),
352  $ work( 2*n+1 ) )
353 *
354  END IF
355 *
356  work( 1 ) = iws
357  RETURN
358 *
359 * End of DGEQP3
360 *
Here is the call graph for this function:
Here is the caller graph for this function:
dlaqp2
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
dnrm2
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dswap
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
dormqr
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:169
dgeqrf
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:147
dlaqps
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