LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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*
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*
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: