LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dlaqps()

subroutine dlaqps ( integer  M,
integer  N,
integer  OFFSET,
integer  NB,
integer  KB,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  JPVT,
double precision, dimension( * )  TAU,
double precision, dimension( * )  VN1,
double precision, dimension( * )  VN2,
double precision, dimension( * )  AUXV,
double precision, dimension( ldf, * )  F,
integer  LDF 
)

DLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3.

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

Purpose:
 DLAQPS computes a step of QR factorization with column pivoting
 of a real M-by-N matrix A by using Blas-3.  It tries to factorize
 NB columns from A starting from the row OFFSET+1, and updates all
 of the matrix with Blas-3 xGEMM.

 In some cases, due to catastrophic cancellations, it cannot
 factorize NB columns.  Hence, the actual number of factorized
 columns is returned in KB.

 Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
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]OFFSET
          OFFSET is INTEGER
          The number of rows of A that have been factorized in
          previous steps.
[in]NB
          NB is INTEGER
          The number of columns to factorize.
[out]KB
          KB is INTEGER
          The number of columns actually factorized.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, block A(OFFSET+1:M,1:KB) is the triangular
          factor obtained and block A(1:OFFSET,1:N) has been
          accordingly pivoted, but no factorized.
          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
          been updated.
[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)
          JPVT(I) = K <==> Column K of the full matrix A has been
          permuted into position I in AP.
[out]TAU
          TAU is DOUBLE PRECISION array, dimension (KB)
          The scalar factors of the elementary reflectors.
[in,out]VN1
          VN1 is DOUBLE PRECISION array, dimension (N)
          The vector with the partial column norms.
[in,out]VN2
          VN2 is DOUBLE PRECISION array, dimension (N)
          The vector with the exact column norms.
[in,out]AUXV
          AUXV is DOUBLE PRECISION array, dimension (NB)
          Auxiliar vector.
[in,out]F
          F is DOUBLE PRECISION array, dimension (LDF,NB)
          Matrix F**T = L*Y**T*A.
[in]LDF
          LDF is INTEGER
          The leading dimension of the array F. LDF >= max(1,N).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors:
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain X. Sun, Computer Science Dept., Duke University, USA
Partial column norm updating strategy modified on April 2011 Z. Drmac and Z. Bujanovic, Dept. of Mathematics, University of Zagreb, Croatia.
References:
LAPACK Working Note 176 [PDF]

Definition at line 179 of file dlaqps.f.

179 *
180 * -- LAPACK auxiliary routine (version 3.7.0) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * December 2016
184 *
185 * .. Scalar Arguments ..
186  INTEGER kb, lda, ldf, m, n, nb, offset
187 * ..
188 * .. Array Arguments ..
189  INTEGER jpvt( * )
190  DOUBLE PRECISION a( lda, * ), auxv( * ), f( ldf, * ), tau( * ),
191  $ vn1( * ), vn2( * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  DOUBLE PRECISION zero, one
198  parameter( zero = 0.0d+0, one = 1.0d+0 )
199 * ..
200 * .. Local Scalars ..
201  INTEGER itemp, j, k, lastrk, lsticc, pvt, rk
202  DOUBLE PRECISION akk, temp, temp2, tol3z
203 * ..
204 * .. External Subroutines ..
205  EXTERNAL dgemm, dgemv, dlarfg, dswap
206 * ..
207 * .. Intrinsic Functions ..
208  INTRINSIC abs, dble, max, min, nint, sqrt
209 * ..
210 * .. External Functions ..
211  INTEGER idamax
212  DOUBLE PRECISION dlamch, dnrm2
213  EXTERNAL idamax, dlamch, dnrm2
214 * ..
215 * .. Executable Statements ..
216 *
217  lastrk = min( m, n+offset )
218  lsticc = 0
219  k = 0
220  tol3z = sqrt(dlamch('Epsilon'))
221 *
222 * Beginning of while loop.
223 *
224  10 CONTINUE
225  IF( ( k.LT.nb ) .AND. ( lsticc.EQ.0 ) ) THEN
226  k = k + 1
227  rk = offset + k
228 *
229 * Determine ith pivot column and swap if necessary
230 *
231  pvt = ( k-1 ) + idamax( n-k+1, vn1( k ), 1 )
232  IF( pvt.NE.k ) THEN
233  CALL dswap( m, a( 1, pvt ), 1, a( 1, k ), 1 )
234  CALL dswap( k-1, f( pvt, 1 ), ldf, f( k, 1 ), ldf )
235  itemp = jpvt( pvt )
236  jpvt( pvt ) = jpvt( k )
237  jpvt( k ) = itemp
238  vn1( pvt ) = vn1( k )
239  vn2( pvt ) = vn2( k )
240  END IF
241 *
242 * Apply previous Householder reflectors to column K:
243 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**T.
244 *
245  IF( k.GT.1 ) THEN
246  CALL dgemv( 'No transpose', m-rk+1, k-1, -one, a( rk, 1 ),
247  $ lda, f( k, 1 ), ldf, one, a( rk, k ), 1 )
248  END IF
249 *
250 * Generate elementary reflector H(k).
251 *
252  IF( rk.LT.m ) THEN
253  CALL dlarfg( m-rk+1, a( rk, k ), a( rk+1, k ), 1, tau( k ) )
254  ELSE
255  CALL dlarfg( 1, a( rk, k ), a( rk, k ), 1, tau( k ) )
256  END IF
257 *
258  akk = a( rk, k )
259  a( rk, k ) = one
260 *
261 * Compute Kth column of F:
262 *
263 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**T*A(RK:M,K).
264 *
265  IF( k.LT.n ) THEN
266  CALL dgemv( 'Transpose', m-rk+1, n-k, tau( k ),
267  $ a( rk, k+1 ), lda, a( rk, k ), 1, zero,
268  $ f( k+1, k ), 1 )
269  END IF
270 *
271 * Padding F(1:K,K) with zeros.
272 *
273  DO 20 j = 1, k
274  f( j, k ) = zero
275  20 CONTINUE
276 *
277 * Incremental updating of F:
278 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**T
279 * *A(RK:M,K).
280 *
281  IF( k.GT.1 ) THEN
282  CALL dgemv( 'Transpose', m-rk+1, k-1, -tau( k ), a( rk, 1 ),
283  $ lda, a( rk, k ), 1, zero, auxv( 1 ), 1 )
284 *
285  CALL dgemv( 'No transpose', n, k-1, one, f( 1, 1 ), ldf,
286  $ auxv( 1 ), 1, one, f( 1, k ), 1 )
287  END IF
288 *
289 * Update the current row of A:
290 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**T.
291 *
292  IF( k.LT.n ) THEN
293  CALL dgemv( 'No transpose', n-k, k, -one, f( k+1, 1 ), ldf,
294  $ a( rk, 1 ), lda, one, a( rk, k+1 ), lda )
295  END IF
296 *
297 * Update partial column norms.
298 *
299  IF( rk.LT.lastrk ) THEN
300  DO 30 j = k + 1, n
301  IF( vn1( j ).NE.zero ) THEN
302 *
303 * NOTE: The following 4 lines follow from the analysis in
304 * Lapack Working Note 176.
305 *
306  temp = abs( a( rk, j ) ) / vn1( j )
307  temp = max( zero, ( one+temp )*( one-temp ) )
308  temp2 = temp*( vn1( j ) / vn2( j ) )**2
309  IF( temp2 .LE. tol3z ) THEN
310  vn2( j ) = dble( lsticc )
311  lsticc = j
312  ELSE
313  vn1( j ) = vn1( j )*sqrt( temp )
314  END IF
315  END IF
316  30 CONTINUE
317  END IF
318 *
319  a( rk, k ) = akk
320 *
321 * End of while loop.
322 *
323  GO TO 10
324  END IF
325  kb = k
326  rk = offset + kb
327 *
328 * Apply the block reflector to the rest of the matrix:
329 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
330 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**T.
331 *
332  IF( kb.LT.min( n, m-offset ) ) THEN
333  CALL dgemm( 'No transpose', 'Transpose', m-rk, n-kb, kb, -one,
334  $ a( rk+1, 1 ), lda, f( kb+1, 1 ), ldf, one,
335  $ a( rk+1, kb+1 ), lda )
336  END IF
337 *
338 * Recomputation of difficult columns.
339 *
340  40 CONTINUE
341  IF( lsticc.GT.0 ) THEN
342  itemp = nint( vn2( lsticc ) )
343  vn1( lsticc ) = dnrm2( m-rk, a( rk+1, lsticc ), 1 )
344 *
345 * NOTE: The computation of VN1( LSTICC ) relies on the fact that
346 * SNRM2 does not fail on vectors with norm below the value of
347 * SQRT(DLAMCH('S'))
348 *
349  vn2( lsticc ) = vn1( lsticc )
350  lsticc = itemp
351  GO TO 40
352  END IF
353 *
354  RETURN
355 *
356 * End of DLAQPS
357 *
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:73
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
Here is the call graph for this function:
Here is the caller graph for this function: