 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.

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).```
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: