LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:106
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition: dznrm2.f90:90
Here is the call graph for this function:
Here is the caller graph for this function: