1 *> \brief \b ZLAQPS computes a step of QR factorization with column pivoting of a real m-by-n matrix A by using BLAS level 3.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZLAQPS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqps.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqps.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqps.f">
21 * SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1,
24 * .. Scalar Arguments ..
25 * INTEGER KB, LDA, LDF, M, N, NB, OFFSET
27 * .. Array Arguments ..
29 * DOUBLE PRECISION VN1( * ), VN2( * )
30 * COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
39 *> ZLAQPS computes a step of QR factorization with column pivoting
40 *> of a complex M-by-N matrix A by using Blas-3. It tries to factorize
41 *> NB columns from A starting from the row OFFSET+1, and updates all
42 *> of the matrix with Blas-3 xGEMM.
44 *> In some cases, due to catastrophic cancellations, it cannot
45 *> factorize NB columns. Hence, the actual number of factorized
46 *> columns is returned in KB.
48 *> Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
57 *> The number of rows of the matrix A. M >= 0.
63 *> The number of columns of the matrix A. N >= 0
69 *> The number of rows of A that have been factorized in
76 *> The number of columns to factorize.
82 *> The number of columns actually factorized.
87 *> A is COMPLEX*16 array, dimension (LDA,N)
88 *> On entry, the M-by-N matrix A.
89 *> On exit, block A(OFFSET+1:M,1:KB) is the triangular
90 *> factor obtained and block A(1:OFFSET,1:N) has been
91 *> accordingly pivoted, but no factorized.
92 *> The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
99 *> The leading dimension of the array A. LDA >= max(1,M).
102 *> \param[in,out] JPVT
104 *> JPVT is INTEGER array, dimension (N)
105 *> JPVT(I) = K <==> Column K of the full matrix A has been
106 *> permuted into position I in AP.
111 *> TAU is COMPLEX*16 array, dimension (KB)
112 *> The scalar factors of the elementary reflectors.
115 *> \param[in,out] VN1
117 *> VN1 is DOUBLE PRECISION array, dimension (N)
118 *> The vector with the partial column norms.
121 *> \param[in,out] VN2
123 *> VN2 is DOUBLE PRECISION array, dimension (N)
124 *> The vector with the exact column norms.
127 *> \param[in,out] AUXV
129 *> AUXV is COMPLEX*16 array, dimension (NB)
135 *> F is COMPLEX*16 array, dimension (LDF,NB)
136 *> Matrix F**H = L * Y**H * A.
142 *> The leading dimension of the array F. LDF >= max(1,N).
148 *> \author Univ. of Tennessee
149 *> \author Univ. of California Berkeley
150 *> \author Univ. of Colorado Denver
153 *> \date September 2012
155 *> \ingroup complex16OTHERauxiliary
157 *> \par Contributors:
160 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
161 *> X. Sun, Computer Science Dept., Duke University, USA
163 *> Partial column norm updating strategy modified on April 2011
164 *> Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
165 *> University of Zagreb, Croatia.
170 *> LAPACK Working Note 176
173 *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a>
176 * =====================================================================
177 SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1,
178 $ VN2, AUXV, F, LDF )
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..--
185 * .. Scalar Arguments ..
186 INTEGER KB, LDA, LDF, M, N, NB, OFFSET
188 * .. Array Arguments ..
190 DOUBLE PRECISION VN1( * ), VN2( * )
191 COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * )
194 * =====================================================================
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 ) )
203 * .. Local Scalars ..
204 INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK
205 DOUBLE PRECISION TEMP, TEMP2, TOL3Z
208 * .. External Subroutines ..
209 EXTERNAL ZGEMM, ZGEMV, ZLARFG, ZSWAP
211 * .. Intrinsic Functions ..
212 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, NINT, SQRT
214 * .. External Functions ..
216 DOUBLE PRECISION DLAMCH, DZNRM2
217 EXTERNAL IDAMAX, DLAMCH, DZNRM2
219 * .. Executable Statements ..
221 LASTRK = MIN( M, N+OFFSET )
224 TOL3Z = SQRT(DLAMCH('Epsilon'))
226 * Beginning of while loop.
229 IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN
233 * Determine ith pivot column and swap if necessary
235 PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 )
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 )
240 JPVT( PVT ) = JPVT( K )
242 VN1( PVT ) = VN1( K )
243 VN2( PVT ) = VN2( K )
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.
251 F( K, J ) = DCONJG( F( K, J ) )
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 )
256 F( K, J ) = DCONJG( F( K, J ) )
260 * Generate elementary reflector H(k).
263 CALL ZLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) )
265 CALL ZLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) )
271 * Compute Kth column of F:
273 * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
276 CALL ZGEMV( 'Conjugate transpose', M-RK+1, N-K, TAU( K ),
277 $ A( RK, K+1 ), LDA, A( RK, K ), 1, CZERO,
281 * Padding F(1:K,K) with zeros.
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
292 CALL ZGEMV( 'Conjugate transpose', M-RK+1, K-1, -TAU( K ),
293 $ A( RK, 1 ), LDA, A( RK, K ), 1, CZERO,
296 CALL ZGEMV( 'No transpose', N, K-1, CONE, F( 1, 1 ), LDF,
297 $ AUXV( 1 ), 1, CONE, F( 1, K ), 1 )
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.
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 )
309 * Update partial column norms.
311 IF( RK.LT.LASTRK ) THEN
313 IF( VN1( J ).NE.ZERO ) THEN
315 * NOTE: The following 4 lines follow from the analysis in
316 * Lapack Working Note 176.
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 )
325 VN1( J ) = VN1( J )*SQRT( TEMP )
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.
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 )
350 * Recomputation of difficult columns.
353 IF( LSTICC.GT.0 ) THEN
354 ITEMP = NINT( VN2( LSTICC ) )
355 VN1( LSTICC ) = DZNRM2( M-RK, A( RK+1, LSTICC ), 1 )
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
361 VN2( LSTICC ) = VN1( LSTICC )