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

◆ zlaqps()

subroutine zlaqps ( integer  m,
integer  n,
integer  offset,
integer  nb,
integer  kb,
complex*16, dimension( lda, * )  a,
integer  lda,
integer, dimension( * )  jpvt,
complex*16, dimension( * )  tau,
double precision, dimension( * )  vn1,
double precision, dimension( * )  vn2,
complex*16, dimension( * )  auxv,
complex*16, dimension( ldf, * )  f,
integer  ldf 
)

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

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

Purpose:
 ZLAQPS 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*16 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*16 array, dimension (KB)
          The scalar factors of the elementary reflectors.
[in,out]VN1
          VN1 is DOUBLE PRECISION array, dimension (N)
          The vector with the partial column norms.
[in,out]VN2
          VN2 is DOUBLE PRECISION array, dimension (N)
          The vector with the exact column norms.
[in,out]AUXV
          AUXV is COMPLEX*16 array, dimension (NB)
          Auxiliary vector.
[in,out]F
          F is COMPLEX*16 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).
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 175 of file zlaqps.f.

177*
178* -- LAPACK auxiliary routine --
179* -- LAPACK is a software package provided by Univ. of Tennessee, --
180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181*
182* .. Scalar Arguments ..
183 INTEGER KB, LDA, LDF, M, N, NB, OFFSET
184* ..
185* .. Array Arguments ..
186 INTEGER JPVT( * )
187 DOUBLE PRECISION VN1( * ), VN2( * )
188 COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
189* ..
190*
191* =====================================================================
192*
193* .. Parameters ..
194 DOUBLE PRECISION ZERO, ONE
195 COMPLEX*16 CZERO, CONE
196 parameter( zero = 0.0d+0, one = 1.0d+0,
197 $ czero = ( 0.0d+0, 0.0d+0 ),
198 $ cone = ( 1.0d+0, 0.0d+0 ) )
199* ..
200* .. Local Scalars ..
201 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
202 DOUBLE PRECISION TEMP, TEMP2, TOL3Z
203 COMPLEX*16 AKK
204* ..
205* .. External Subroutines ..
206 EXTERNAL zgemm, zgemv, zlarfg, zswap
207* ..
208* .. Intrinsic Functions ..
209 INTRINSIC abs, dble, dconjg, max, min, nint, sqrt
210* ..
211* .. External Functions ..
212 INTEGER IDAMAX
213 DOUBLE PRECISION DLAMCH, DZNRM2
214 EXTERNAL idamax, dlamch, dznrm2
215* ..
216* .. Executable Statements ..
217*
218 lastrk = min( m, n+offset )
219 lsticc = 0
220 k = 0
221 tol3z = sqrt(dlamch('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 ) + idamax( n-k+1, vn1( k ), 1 )
233 IF( pvt.NE.k ) THEN
234 CALL zswap( m, a( 1, pvt ), 1, a( 1, k ), 1 )
235 CALL zswap( 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)**H.
245*
246 IF( k.GT.1 ) THEN
247 DO 20 j = 1, k - 1
248 f( k, j ) = dconjg( f( k, j ) )
249 20 CONTINUE
250 CALL zgemv( 'No transpose', m-rk+1, k-1, -cone, a( rk, 1 ),
251 $ lda, f( k, 1 ), ldf, cone, a( rk, k ), 1 )
252 DO 30 j = 1, k - 1
253 f( k, j ) = dconjg( f( k, j ) )
254 30 CONTINUE
255 END IF
256*
257* Generate elementary reflector H(k).
258*
259 IF( rk.LT.m ) THEN
260 CALL zlarfg( m-rk+1, a( rk, k ), a( rk+1, k ), 1, tau( k ) )
261 ELSE
262 CALL zlarfg( 1, a( rk, k ), a( rk, k ), 1, tau( k ) )
263 END IF
264*
265 akk = a( rk, k )
266 a( rk, k ) = cone
267*
268* Compute Kth column of F:
269*
270* Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
271*
272 IF( k.LT.n ) THEN
273 CALL zgemv( 'Conjugate transpose', m-rk+1, n-k, tau( k ),
274 $ a( rk, k+1 ), lda, a( rk, k ), 1, czero,
275 $ f( k+1, k ), 1 )
276 END IF
277*
278* Padding F(1:K,K) with zeros.
279*
280 DO 40 j = 1, k
281 f( j, k ) = czero
282 40 CONTINUE
283*
284* Incremental updating of F:
285* F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H
286* *A(RK:M,K).
287*
288 IF( k.GT.1 ) THEN
289 CALL zgemv( 'Conjugate transpose', m-rk+1, k-1, -tau( k ),
290 $ a( rk, 1 ), lda, a( rk, k ), 1, czero,
291 $ auxv( 1 ), 1 )
292*
293 CALL zgemv( 'No transpose', n, k-1, cone, f( 1, 1 ), ldf,
294 $ auxv( 1 ), 1, cone, f( 1, k ), 1 )
295 END IF
296*
297* Update the current row of A:
298* A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**H.
299*
300 IF( k.LT.n ) THEN
301 CALL zgemm( 'No transpose', 'Conjugate transpose', 1, n-k,
302 $ k, -cone, a( rk, 1 ), lda, f( k+1, 1 ), ldf,
303 $ cone, a( rk, k+1 ), lda )
304 END IF
305*
306* Update partial column norms.
307*
308 IF( rk.LT.lastrk ) THEN
309 DO 50 j = k + 1, n
310 IF( vn1( j ).NE.zero ) THEN
311*
312* NOTE: The following 4 lines follow from the analysis in
313* Lapack Working Note 176.
314*
315 temp = abs( a( rk, j ) ) / vn1( j )
316 temp = max( zero, ( one+temp )*( one-temp ) )
317 temp2 = temp*( vn1( j ) / vn2( j ) )**2
318 IF( temp2 .LE. tol3z ) THEN
319 vn2( j ) = dble( lsticc )
320 lsticc = j
321 ELSE
322 vn1( j ) = vn1( j )*sqrt( temp )
323 END IF
324 END IF
325 50 CONTINUE
326 END IF
327*
328 a( rk, k ) = akk
329*
330* End of while loop.
331*
332 GO TO 10
333 END IF
334 kb = k
335 rk = offset + kb
336*
337* Apply the block reflector to the rest of the matrix:
338* A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
339* A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**H.
340*
341 IF( kb.LT.min( n, m-offset ) ) THEN
342 CALL zgemm( 'No transpose', 'Conjugate transpose', m-rk, n-kb,
343 $ kb, -cone, a( rk+1, 1 ), lda, f( kb+1, 1 ), ldf,
344 $ cone, a( rk+1, kb+1 ), lda )
345 END IF
346*
347* Recomputation of difficult columns.
348*
349 60 CONTINUE
350 IF( lsticc.GT.0 ) THEN
351 itemp = nint( vn2( lsticc ) )
352 vn1( lsticc ) = dznrm2( m-rk, a( rk+1, lsticc ), 1 )
353*
354* NOTE: The computation of VN1( LSTICC ) relies on the fact that
355* SNRM2 does not fail on vectors with norm below the value of
356* SQRT(DLAMCH('S'))
357*
358 vn2( lsticc ) = vn1( lsticc )
359 lsticc = itemp
360 GO TO 60
361 END IF
362*
363 RETURN
364*
365* End of ZLAQPS
366*
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlarfg
subroutine zlarfg(n, alpha, x, incx, tau)
ZLARFG generates an elementary reflector (Householder matrix).
Definition zlarfg.f:106
dznrm2
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: