LAPACK  3.10.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)
          Auxiliary 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.
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 176 of file slaqps.f.

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