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