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. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLAQPS + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqps.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqps.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqps.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, 22* VN2, AUXV, F, LDF ) 23* 24* .. Scalar Arguments .. 25* INTEGER KB, LDA, LDF, M, N, NB, OFFSET 26* .. 27* .. Array Arguments .. 28* INTEGER JPVT( * ) 29* DOUBLE PRECISION VN1( * ), VN2( * ) 30* COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 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. 43*> 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. 47*> 48*> Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. 49*> \endverbatim 50* 51* Arguments: 52* ========== 53* 54*> \param[in] M 55*> \verbatim 56*> M is INTEGER 57*> The number of rows of the matrix A. M >= 0. 58*> \endverbatim 59*> 60*> \param[in] N 61*> \verbatim 62*> N is INTEGER 63*> The number of columns of the matrix A. N >= 0 64*> \endverbatim 65*> 66*> \param[in] OFFSET 67*> \verbatim 68*> OFFSET is INTEGER 69*> The number of rows of A that have been factorized in 70*> previous steps. 71*> \endverbatim 72*> 73*> \param[in] NB 74*> \verbatim 75*> NB is INTEGER 76*> The number of columns to factorize. 77*> \endverbatim 78*> 79*> \param[out] KB 80*> \verbatim 81*> KB is INTEGER 82*> The number of columns actually factorized. 83*> \endverbatim 84*> 85*> \param[in,out] A 86*> \verbatim 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 93*> been updated. 94*> \endverbatim 95*> 96*> \param[in] LDA 97*> \verbatim 98*> LDA is INTEGER 99*> The leading dimension of the array A. LDA >= max(1,M). 100*> \endverbatim 101*> 102*> \param[in,out] JPVT 103*> \verbatim 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. 107*> \endverbatim 108*> 109*> \param[out] TAU 110*> \verbatim 111*> TAU is COMPLEX*16 array, dimension (KB) 112*> The scalar factors of the elementary reflectors. 113*> \endverbatim 114*> 115*> \param[in,out] VN1 116*> \verbatim 117*> VN1 is DOUBLE PRECISION array, dimension (N) 118*> The vector with the partial column norms. 119*> \endverbatim 120*> 121*> \param[in,out] VN2 122*> \verbatim 123*> VN2 is DOUBLE PRECISION array, dimension (N) 124*> The vector with the exact column norms. 125*> \endverbatim 126*> 127*> \param[in,out] AUXV 128*> \verbatim 129*> AUXV is COMPLEX*16 array, dimension (NB) 130*> Auxiliary vector. 131*> \endverbatim 132*> 133*> \param[in,out] F 134*> \verbatim 135*> F is COMPLEX*16 array, dimension (LDF,NB) 136*> Matrix F**H = L * Y**H * A. 137*> \endverbatim 138*> 139*> \param[in] LDF 140*> \verbatim 141*> LDF is INTEGER 142*> The leading dimension of the array F. LDF >= max(1,N). 143*> \endverbatim 144* 145* Authors: 146* ======== 147* 148*> \author Univ. of Tennessee 149*> \author Univ. of California Berkeley 150*> \author Univ. of Colorado Denver 151*> \author NAG Ltd. 152* 153*> \ingroup complex16OTHERauxiliary 154* 155*> \par Contributors: 156* ================== 157*> 158*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 159*> X. Sun, Computer Science Dept., Duke University, USA 160*> \n 161*> Partial column norm updating strategy modified on April 2011 162*> Z. Drmac and Z. Bujanovic, Dept. of Mathematics, 163*> University of Zagreb, Croatia. 164* 165*> \par References: 166* ================ 167*> 168*> LAPACK Working Note 176 169* 170*> \htmlonly 171*> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a> 172*> \endhtmlonly 173* 174* ===================================================================== 175 SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, 176 $ VN2, AUXV, F, LDF ) 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* 367 END 368