LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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)
          Auxiliar 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.
Date
September 2012
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 179 of file zlaqps.f.

179 *
180 * -- LAPACK auxiliary routine (version 3.4.2) --
181 * -- LAPACK is a software package provided by Univ. of Tennessee, --
182 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * September 2012
184 *
185 * .. Scalar Arguments ..
186  INTEGER kb, lda, ldf, m, n, nb, offset
187 * ..
188 * .. Array Arguments ..
189  INTEGER jpvt( * )
190  DOUBLE PRECISION vn1( * ), vn2( * )
191  COMPLEX*16 a( lda, * ), auxv( * ), f( ldf, * ), tau( * )
192 * ..
193 *
194 * =====================================================================
195 *
196 * .. Parameters ..
197  DOUBLE PRECISION zero, one
198  COMPLEX*16 czero, cone
199  parameter ( zero = 0.0d+0, one = 1.0d+0,
200  $ czero = ( 0.0d+0, 0.0d+0 ),
201  $ cone = ( 1.0d+0, 0.0d+0 ) )
202 * ..
203 * .. Local Scalars ..
204  INTEGER itemp, j, k, lastrk, lsticc, pvt, rk
205  DOUBLE PRECISION temp, temp2, tol3z
206  COMPLEX*16 akk
207 * ..
208 * .. External Subroutines ..
209  EXTERNAL zgemm, zgemv, zlarfg, zswap
210 * ..
211 * .. Intrinsic Functions ..
212  INTRINSIC abs, dble, dconjg, max, min, nint, sqrt
213 * ..
214 * .. External Functions ..
215  INTEGER idamax
216  DOUBLE PRECISION dlamch, dznrm2
217  EXTERNAL idamax, dlamch, dznrm2
218 * ..
219 * .. Executable Statements ..
220 *
221  lastrk = min( m, n+offset )
222  lsticc = 0
223  k = 0
224  tol3z = sqrt(dlamch('Epsilon'))
225 *
226 * Beginning of while loop.
227 *
228  10 CONTINUE
229  IF( ( k.LT.nb ) .AND. ( lsticc.EQ.0 ) ) THEN
230  k = k + 1
231  rk = offset + k
232 *
233 * Determine ith pivot column and swap if necessary
234 *
235  pvt = ( k-1 ) + idamax( n-k+1, vn1( k ), 1 )
236  IF( pvt.NE.k ) THEN
237  CALL zswap( m, a( 1, pvt ), 1, a( 1, k ), 1 )
238  CALL zswap( k-1, f( pvt, 1 ), ldf, f( k, 1 ), ldf )
239  itemp = jpvt( pvt )
240  jpvt( pvt ) = jpvt( k )
241  jpvt( k ) = itemp
242  vn1( pvt ) = vn1( k )
243  vn2( pvt ) = vn2( k )
244  END IF
245 *
246 * Apply previous Householder reflectors to column K:
247 * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**H.
248 *
249  IF( k.GT.1 ) THEN
250  DO 20 j = 1, k - 1
251  f( k, j ) = dconjg( f( k, j ) )
252  20 CONTINUE
253  CALL zgemv( 'No transpose', m-rk+1, k-1, -cone, a( rk, 1 ),
254  $ lda, f( k, 1 ), ldf, cone, a( rk, k ), 1 )
255  DO 30 j = 1, k - 1
256  f( k, j ) = dconjg( f( k, j ) )
257  30 CONTINUE
258  END IF
259 *
260 * Generate elementary reflector H(k).
261 *
262  IF( rk.LT.m ) THEN
263  CALL zlarfg( m-rk+1, a( rk, k ), a( rk+1, k ), 1, tau( k ) )
264  ELSE
265  CALL zlarfg( 1, a( rk, k ), a( rk, k ), 1, tau( k ) )
266  END IF
267 *
268  akk = a( rk, k )
269  a( rk, k ) = cone
270 *
271 * Compute Kth column of F:
272 *
273 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
274 *
275  IF( k.LT.n ) THEN
276  CALL zgemv( 'Conjugate transpose', m-rk+1, n-k, tau( k ),
277  $ a( rk, k+1 ), lda, a( rk, k ), 1, czero,
278  $ f( k+1, k ), 1 )
279  END IF
280 *
281 * Padding F(1:K,K) with zeros.
282 *
283  DO 40 j = 1, k
284  f( j, k ) = czero
285  40 CONTINUE
286 *
287 * Incremental updating of F:
288 * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H
289 * *A(RK:M,K).
290 *
291  IF( k.GT.1 ) THEN
292  CALL zgemv( 'Conjugate transpose', m-rk+1, k-1, -tau( k ),
293  $ a( rk, 1 ), lda, a( rk, k ), 1, czero,
294  $ auxv( 1 ), 1 )
295 *
296  CALL zgemv( 'No transpose', n, k-1, cone, f( 1, 1 ), ldf,
297  $ auxv( 1 ), 1, cone, f( 1, k ), 1 )
298  END IF
299 *
300 * Update the current row of A:
301 * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**H.
302 *
303  IF( k.LT.n ) THEN
304  CALL zgemm( 'No transpose', 'Conjugate transpose', 1, n-k,
305  $ k, -cone, a( rk, 1 ), lda, f( k+1, 1 ), ldf,
306  $ cone, a( rk, k+1 ), lda )
307  END IF
308 *
309 * Update partial column norms.
310 *
311  IF( rk.LT.lastrk ) THEN
312  DO 50 j = k + 1, n
313  IF( vn1( j ).NE.zero ) THEN
314 *
315 * NOTE: The following 4 lines follow from the analysis in
316 * Lapack Working Note 176.
317 *
318  temp = abs( a( rk, j ) ) / vn1( j )
319  temp = max( zero, ( one+temp )*( one-temp ) )
320  temp2 = temp*( vn1( j ) / vn2( j ) )**2
321  IF( temp2 .LE. tol3z ) THEN
322  vn2( j ) = dble( lsticc )
323  lsticc = j
324  ELSE
325  vn1( j ) = vn1( j )*sqrt( temp )
326  END IF
327  END IF
328  50 CONTINUE
329  END IF
330 *
331  a( rk, k ) = akk
332 *
333 * End of while loop.
334 *
335  GO TO 10
336  END IF
337  kb = k
338  rk = offset + kb
339 *
340 * Apply the block reflector to the rest of the matrix:
341 * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
342 * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**H.
343 *
344  IF( kb.LT.min( n, m-offset ) ) THEN
345  CALL zgemm( 'No transpose', 'Conjugate transpose', m-rk, n-kb,
346  $ kb, -cone, a( rk+1, 1 ), lda, f( kb+1, 1 ), ldf,
347  $ cone, a( rk+1, kb+1 ), lda )
348  END IF
349 *
350 * Recomputation of difficult columns.
351 *
352  60 CONTINUE
353  IF( lsticc.GT.0 ) THEN
354  itemp = nint( vn2( lsticc ) )
355  vn1( lsticc ) = dznrm2( m-rk, a( rk+1, lsticc ), 1 )
356 *
357 * NOTE: The computation of VN1( LSTICC ) relies on the fact that
358 * SNRM2 does not fail on vectors with norm below the value of
359 * SQRT(DLAMCH('S'))
360 *
361  vn2( lsticc ) = vn1( lsticc )
362  lsticc = itemp
363  GO TO 60
364  END IF
365 *
366  RETURN
367 *
368 * End of ZLAQPS
369 *
idamax
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zlarfg
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
Definition: zlarfg.f:108
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
dznrm2
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56

Here is the call graph for this function:

Here is the caller graph for this function: