LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

## ◆ claqps()

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

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

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

Purpose:
``` CLAQPS computes a step of QR factorization with column pivoting
of a complex 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 COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (NB) Auxiliary vector.``` [in,out] F ``` F is COMPLEX array, dimension (LDF,NB) Matrix F**H = L * Y**H * 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 claqps.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 VN1( * ), VN2( * )
189 COMPLEX A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
190* ..
191*
192* =====================================================================
193*
194* .. Parameters ..
195 REAL ZERO, ONE
196 COMPLEX CZERO, CONE
197 parameter( zero = 0.0e+0, one = 1.0e+0,
198 \$ czero = ( 0.0e+0, 0.0e+0 ),
199 \$ cone = ( 1.0e+0, 0.0e+0 ) )
200* ..
201* .. Local Scalars ..
202 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
203 REAL TEMP, TEMP2, TOL3Z
204 COMPLEX AKK
205* ..
206* .. External Subroutines ..
207 EXTERNAL cgemm, cgemv, clarfg, cswap
208* ..
209* .. Intrinsic Functions ..
210 INTRINSIC abs, conjg, max, min, nint, real, sqrt
211* ..
212* .. External Functions ..
213 INTEGER ISAMAX
214 REAL SCNRM2, SLAMCH
215 EXTERNAL isamax, scnrm2, slamch
216* ..
217* .. Executable Statements ..
218*
219 lastrk = min( m, n+offset )
220 lsticc = 0
221 k = 0
222 tol3z = sqrt(slamch('Epsilon'))
223*
224* Beginning of while loop.
225*
226 10 CONTINUE
227 IF( ( k.LT.nb ) .AND. ( lsticc.EQ.0 ) ) THEN
228 k = k + 1
229 rk = offset + k
230*
231* Determine ith pivot column and swap if necessary
232*
233 pvt = ( k-1 ) + isamax( n-k+1, vn1( k ), 1 )
234 IF( pvt.NE.k ) THEN
235 CALL cswap( m, a( 1, pvt ), 1, a( 1, k ), 1 )
236 CALL cswap( k-1, f( pvt, 1 ), ldf, f( k, 1 ), ldf )
237 itemp = jpvt( pvt )
238 jpvt( pvt ) = jpvt( k )
239 jpvt( k ) = itemp
240 vn1( pvt ) = vn1( k )
241 vn2( pvt ) = vn2( k )
242 END IF
243*
244* Apply previous Householder reflectors to column K:
245* A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**H.
246*
247 IF( k.GT.1 ) THEN
248 DO 20 j = 1, k - 1
249 f( k, j ) = conjg( f( k, j ) )
250 20 CONTINUE
251 CALL cgemv( 'No transpose', m-rk+1, k-1, -cone, a( rk, 1 ),
252 \$ lda, f( k, 1 ), ldf, cone, a( rk, k ), 1 )
253 DO 30 j = 1, k - 1
254 f( k, j ) = conjg( f( k, j ) )
255 30 CONTINUE
256 END IF
257*
258* Generate elementary reflector H(k).
259*
260 IF( rk.LT.m ) THEN
261 CALL clarfg( m-rk+1, a( rk, k ), a( rk+1, k ), 1, tau( k ) )
262 ELSE
263 CALL clarfg( 1, a( rk, k ), a( rk, k ), 1, tau( k ) )
264 END IF
265*
266 akk = a( rk, k )
267 a( rk, k ) = cone
268*
269* Compute Kth column of F:
270*
271* Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
272*
273 IF( k.LT.n ) THEN
274 CALL cgemv( 'Conjugate transpose', m-rk+1, n-k, tau( k ),
275 \$ a( rk, k+1 ), lda, a( rk, k ), 1, czero,
276 \$ f( k+1, k ), 1 )
277 END IF
278*
279* Padding F(1:K,K) with zeros.
280*
281 DO 40 j = 1, k
282 f( j, k ) = czero
283 40 CONTINUE
284*
285* Incremental updating of F:
286* F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H
287* *A(RK:M,K).
288*
289 IF( k.GT.1 ) THEN
290 CALL cgemv( 'Conjugate transpose', m-rk+1, k-1, -tau( k ),
291 \$ a( rk, 1 ), lda, a( rk, k ), 1, czero,
292 \$ auxv( 1 ), 1 )
293*
294 CALL cgemv( 'No transpose', n, k-1, cone, f( 1, 1 ), ldf,
295 \$ auxv( 1 ), 1, cone, f( 1, k ), 1 )
296 END IF
297*
298* Update the current row of A:
299* A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**H.
300*
301 IF( k.LT.n ) THEN
302 CALL cgemm( 'No transpose', 'Conjugate transpose', 1, n-k,
303 \$ k, -cone, a( rk, 1 ), lda, f( k+1, 1 ), ldf,
304 \$ cone, a( rk, k+1 ), lda )
305 END IF
306*
307* Update partial column norms.
308*
309 IF( rk.LT.lastrk ) THEN
310 DO 50 j = k + 1, n
311 IF( vn1( j ).NE.zero ) THEN
312*
313* NOTE: The following 4 lines follow from the analysis in
314* Lapack Working Note 176.
315*
316 temp = abs( a( rk, j ) ) / vn1( j )
317 temp = max( zero, ( one+temp )*( one-temp ) )
318 temp2 = temp*( vn1( j ) / vn2( j ) )**2
319 IF( temp2 .LE. tol3z ) THEN
320 vn2( j ) = real( lsticc )
321 lsticc = j
322 ELSE
323 vn1( j ) = vn1( j )*sqrt( temp )
324 END IF
325 END IF
326 50 CONTINUE
327 END IF
328*
329 a( rk, k ) = akk
330*
331* End of while loop.
332*
333 GO TO 10
334 END IF
335 kb = k
336 rk = offset + kb
337*
338* Apply the block reflector to the rest of the matrix:
339* A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
340* A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**H.
341*
342 IF( kb.LT.min( n, m-offset ) ) THEN
343 CALL cgemm( 'No transpose', 'Conjugate transpose', m-rk, n-kb,
344 \$ kb, -cone, a( rk+1, 1 ), lda, f( kb+1, 1 ), ldf,
345 \$ cone, a( rk+1, kb+1 ), lda )
346 END IF
347*
348* Recomputation of difficult columns.
349*
350 60 CONTINUE
351 IF( lsticc.GT.0 ) THEN
352 itemp = nint( vn2( lsticc ) )
353 vn1( lsticc ) = scnrm2( m-rk, a( rk+1, lsticc ), 1 )
354*
355* NOTE: The computation of VN1( LSTICC ) relies on the fact that
356* SNRM2 does not fail on vectors with norm below the value of
357* SQRT(DLAMCH('S'))
358*
359 vn2( lsticc ) = vn1( lsticc )
360 lsticc = itemp
361 GO TO 60
362 END IF
363*
364 RETURN
365*
366* End of CLAQPS
367*
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:106
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
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: