LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ slaqps()

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

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

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

Purpose:
 SLAQPS 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 REAL 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 REAL array, dimension (KB)
          The scalar factors of the elementary reflectors.
[in,out]VN1
          VN1 is REAL array, dimension (N)
          The vector with the partial column norms.
[in,out]VN2
          VN2 is REAL array, dimension (N)
          The vector with the exact column norms.
[in,out]AUXV
          AUXV is REAL array, dimension (NB)
          Auxiliar vector.
[in,out]F
          F is REAL 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 180 of file slaqps.f.

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