LAPACK 3.12.1
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 171 of file zlaqps.f.

174*
175* -- LAPACK auxiliary routine --
176* -- LAPACK is a software package provided by Univ. of Tennessee, --
177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*
179* .. Scalar Arguments ..
180 INTEGER KB, LDA, LDF, M, N, NB, OFFSET
181* ..
182* .. Array Arguments ..
183 INTEGER JPVT( * )
184 DOUBLE PRECISION VN1( * ), VN2( * )
185 COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ZERO, ONE
192 COMPLEX*16 CZERO, CONE
193 parameter( zero = 0.0d+0, one = 1.0d+0,
194 $ czero = ( 0.0d+0, 0.0d+0 ),
195 $ cone = ( 1.0d+0, 0.0d+0 ) )
196* ..
197* .. Local Scalars ..
198 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
199 DOUBLE PRECISION TEMP, TEMP2, TOL3Z
200 COMPLEX*16 AKK
201* ..
202* .. External Subroutines ..
203 EXTERNAL zgemm, zgemv, zlarfg, zswap
204* ..
205* .. Intrinsic Functions ..
206 INTRINSIC abs, dble, dconjg, max, min, nint, sqrt
207* ..
208* .. External Functions ..
209 INTEGER IDAMAX
210 DOUBLE PRECISION DLAMCH, DZNRM2
211 EXTERNAL idamax, dlamch, dznrm2
212* ..
213* .. Executable Statements ..
214*
215 lastrk = min( m, n+offset )
216 lsticc = 0
217 k = 0
218 tol3z = sqrt(dlamch('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 ) + idamax( n-k+1, vn1( k ), 1 )
230 IF( pvt.NE.k ) THEN
231 CALL zswap( m, a( 1, pvt ), 1, a( 1, k ), 1 )
232 CALL zswap( 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)**H.
242*
243 IF( k.GT.1 ) THEN
244 DO 20 j = 1, k - 1
245 f( k, j ) = dconjg( f( k, j ) )
246 20 CONTINUE
247 CALL zgemv( 'No transpose', m-rk+1, k-1, -cone, a( rk,
248 $ 1 ),
249 $ lda, f( k, 1 ), ldf, cone, a( rk, k ), 1 )
250 DO 30 j = 1, k - 1
251 f( k, j ) = dconjg( f( k, j ) )
252 30 CONTINUE
253 END IF
254*
255* Generate elementary reflector H(k).
256*
257 IF( rk.LT.m ) THEN
258 CALL zlarfg( m-rk+1, a( rk, k ), a( rk+1, k ), 1,
259 $ tau( k ) )
260 ELSE
261 CALL zlarfg( 1, a( rk, k ), a( rk, k ), 1, tau( k ) )
262 END IF
263*
264 akk = a( rk, k )
265 a( rk, k ) = cone
266*
267* Compute Kth column of F:
268*
269* Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
270*
271 IF( k.LT.n ) THEN
272 CALL zgemv( 'Conjugate transpose', m-rk+1, n-k, tau( k ),
273 $ a( rk, k+1 ), lda, a( rk, k ), 1, czero,
274 $ f( k+1, k ), 1 )
275 END IF
276*
277* Padding F(1:K,K) with zeros.
278*
279 DO 40 j = 1, k
280 f( j, k ) = czero
281 40 CONTINUE
282*
283* Incremental updating of F:
284* F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H
285* *A(RK:M,K).
286*
287 IF( k.GT.1 ) THEN
288 CALL zgemv( 'Conjugate transpose', m-rk+1, k-1,
289 $ -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,
302 $ 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 ) = dble( 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 zgemm( 'No transpose', 'Conjugate transpose', m-rk,
344 $ n-kb,
345 $ kb, -cone, a( rk+1, 1 ), lda, f( kb+1, 1 ), ldf,
346 $ cone, a( rk+1, kb+1 ), lda )
347 END IF
348*
349* Recomputation of difficult columns.
350*
351 60 CONTINUE
352 IF( lsticc.GT.0 ) THEN
353 itemp = nint( vn2( lsticc ) )
354 vn1( lsticc ) = dznrm2( m-rk, a( rk+1, lsticc ), 1 )
355*
356* NOTE: The computation of VN1( LSTICC ) relies on the fact that
357* SNRM2 does not fail on vectors with norm below the value of
358* SQRT(DLAMCH('S'))
359*
360 vn2( lsticc ) = vn1( lsticc )
361 lsticc = itemp
362 GO TO 60
363 END IF
364*
365 RETURN
366*
367* End of ZLAQPS
368*
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:104
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: