LAPACK  3.10.1
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)
          Auxiliary 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.
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 175 of file dlaqps.f.

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