![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, 2: $ VN2, AUXV, F, LDF ) 3: * 4: * -- LAPACK auxiliary routine (version 3.2.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * June 2010 8: * 9: * .. Scalar Arguments .. 10: INTEGER KB, LDA, LDF, M, N, NB, OFFSET 11: * .. 12: * .. Array Arguments .. 13: INTEGER JPVT( * ) 14: DOUBLE PRECISION VN1( * ), VN2( * ) 15: COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * ZLAQPS computes a step of QR factorization with column pivoting 22: * of a complex M-by-N matrix A by using Blas-3. It tries to factorize 23: * NB columns from A starting from the row OFFSET+1, and updates all 24: * of the matrix with Blas-3 xGEMM. 25: * 26: * In some cases, due to catastrophic cancellations, it cannot 27: * factorize NB columns. Hence, the actual number of factorized 28: * columns is returned in KB. 29: * 30: * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. 31: * 32: * Arguments 33: * ========= 34: * 35: * M (input) INTEGER 36: * The number of rows of the matrix A. M >= 0. 37: * 38: * N (input) INTEGER 39: * The number of columns of the matrix A. N >= 0 40: * 41: * OFFSET (input) INTEGER 42: * The number of rows of A that have been factorized in 43: * previous steps. 44: * 45: * NB (input) INTEGER 46: * The number of columns to factorize. 47: * 48: * KB (output) INTEGER 49: * The number of columns actually factorized. 50: * 51: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 52: * On entry, the M-by-N matrix A. 53: * On exit, block A(OFFSET+1:M,1:KB) is the triangular 54: * factor obtained and block A(1:OFFSET,1:N) has been 55: * accordingly pivoted, but no factorized. 56: * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has 57: * been updated. 58: * 59: * LDA (input) INTEGER 60: * The leading dimension of the array A. LDA >= max(1,M). 61: * 62: * JPVT (input/output) INTEGER array, dimension (N) 63: * JPVT(I) = K <==> Column K of the full matrix A has been 64: * permuted into position I in AP. 65: * 66: * TAU (output) COMPLEX*16 array, dimension (KB) 67: * The scalar factors of the elementary reflectors. 68: * 69: * VN1 (input/output) DOUBLE PRECISION array, dimension (N) 70: * The vector with the partial column norms. 71: * 72: * VN2 (input/output) DOUBLE PRECISION array, dimension (N) 73: * The vector with the exact column norms. 74: * 75: * AUXV (input/output) COMPLEX*16 array, dimension (NB) 76: * Auxiliar vector. 77: * 78: * F (input/output) COMPLEX*16 array, dimension (LDF,NB) 79: * Matrix F' = L*Y'*A. 80: * 81: * LDF (input) INTEGER 82: * The leading dimension of the array F. LDF >= max(1,N). 83: * 84: * Further Details 85: * =============== 86: * 87: * Based on contributions by 88: * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 89: * X. Sun, Computer Science Dept., Duke University, USA 90: * 91: * ===================================================================== 92: * 93: * .. Parameters .. 94: DOUBLE PRECISION ZERO, ONE 95: COMPLEX*16 CZERO, CONE 96: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 97: $ CZERO = ( 0.0D+0, 0.0D+0 ), 98: $ CONE = ( 1.0D+0, 0.0D+0 ) ) 99: * .. 100: * .. Local Scalars .. 101: INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK 102: DOUBLE PRECISION TEMP, TEMP2, TOL3Z 103: COMPLEX*16 AKK 104: * .. 105: * .. External Subroutines .. 106: EXTERNAL ZGEMM, ZGEMV, ZLARFG, ZSWAP 107: * .. 108: * .. Intrinsic Functions .. 109: INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, NINT, SQRT 110: * .. 111: * .. External Functions .. 112: INTEGER IDAMAX 113: DOUBLE PRECISION DLAMCH, DZNRM2 114: EXTERNAL IDAMAX, DLAMCH, DZNRM2 115: * .. 116: * .. Executable Statements .. 117: * 118: LASTRK = MIN( M, N+OFFSET ) 119: LSTICC = 0 120: K = 0 121: TOL3Z = SQRT(DLAMCH('Epsilon')) 122: * 123: * Beginning of while loop. 124: * 125: 10 CONTINUE 126: IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN 127: K = K + 1 128: RK = OFFSET + K 129: * 130: * Determine ith pivot column and swap if necessary 131: * 132: PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 ) 133: IF( PVT.NE.K ) THEN 134: CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) 135: CALL ZSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) 136: ITEMP = JPVT( PVT ) 137: JPVT( PVT ) = JPVT( K ) 138: JPVT( K ) = ITEMP 139: VN1( PVT ) = VN1( K ) 140: VN2( PVT ) = VN2( K ) 141: END IF 142: * 143: * Apply previous Householder reflectors to column K: 144: * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'. 145: * 146: IF( K.GT.1 ) THEN 147: DO 20 J = 1, K - 1 148: F( K, J ) = DCONJG( F( K, J ) ) 149: 20 CONTINUE 150: CALL ZGEMV( 'No transpose', M-RK+1, K-1, -CONE, A( RK, 1 ), 151: $ LDA, F( K, 1 ), LDF, CONE, A( RK, K ), 1 ) 152: DO 30 J = 1, K - 1 153: F( K, J ) = DCONJG( F( K, J ) ) 154: 30 CONTINUE 155: END IF 156: * 157: * Generate elementary reflector H(k). 158: * 159: IF( RK.LT.M ) THEN 160: CALL ZLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) 161: ELSE 162: CALL ZLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) 163: END IF 164: * 165: AKK = A( RK, K ) 166: A( RK, K ) = CONE 167: * 168: * Compute Kth column of F: 169: * 170: * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K). 171: * 172: IF( K.LT.N ) THEN 173: CALL ZGEMV( 'Conjugate transpose', M-RK+1, N-K, TAU( K ), 174: $ A( RK, K+1 ), LDA, A( RK, K ), 1, CZERO, 175: $ F( K+1, K ), 1 ) 176: END IF 177: * 178: * Padding F(1:K,K) with zeros. 179: * 180: DO 40 J = 1, K 181: F( J, K ) = CZERO 182: 40 CONTINUE 183: * 184: * Incremental updating of F: 185: * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)' 186: * *A(RK:M,K). 187: * 188: IF( K.GT.1 ) THEN 189: CALL ZGEMV( 'Conjugate transpose', M-RK+1, K-1, -TAU( K ), 190: $ A( RK, 1 ), LDA, A( RK, K ), 1, CZERO, 191: $ AUXV( 1 ), 1 ) 192: * 193: CALL ZGEMV( 'No transpose', N, K-1, CONE, F( 1, 1 ), LDF, 194: $ AUXV( 1 ), 1, CONE, F( 1, K ), 1 ) 195: END IF 196: * 197: * Update the current row of A: 198: * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'. 199: * 200: IF( K.LT.N ) THEN 201: CALL ZGEMM( 'No transpose', 'Conjugate transpose', 1, N-K, 202: $ K, -CONE, A( RK, 1 ), LDA, F( K+1, 1 ), LDF, 203: $ CONE, A( RK, K+1 ), LDA ) 204: END IF 205: * 206: * Update partial column norms. 207: * 208: IF( RK.LT.LASTRK ) THEN 209: DO 50 J = K + 1, N 210: IF( VN1( J ).NE.ZERO ) THEN 211: * 212: * NOTE: The following 4 lines follow from the analysis in 213: * Lapack Working Note 176. 214: * 215: TEMP = ABS( A( RK, J ) ) / VN1( J ) 216: TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) 217: TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 218: IF( TEMP2 .LE. TOL3Z ) THEN 219: VN2( J ) = DBLE( LSTICC ) 220: LSTICC = J 221: ELSE 222: VN1( J ) = VN1( J )*SQRT( TEMP ) 223: END IF 224: END IF 225: 50 CONTINUE 226: END IF 227: * 228: A( RK, K ) = AKK 229: * 230: * End of while loop. 231: * 232: GO TO 10 233: END IF 234: KB = K 235: RK = OFFSET + KB 236: * 237: * Apply the block reflector to the rest of the matrix: 238: * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - 239: * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'. 240: * 241: IF( KB.LT.MIN( N, M-OFFSET ) ) THEN 242: CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-RK, N-KB, 243: $ KB, -CONE, A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, 244: $ CONE, A( RK+1, KB+1 ), LDA ) 245: END IF 246: * 247: * Recomputation of difficult columns. 248: * 249: 60 CONTINUE 250: IF( LSTICC.GT.0 ) THEN 251: ITEMP = NINT( VN2( LSTICC ) ) 252: VN1( LSTICC ) = DZNRM2( M-RK, A( RK+1, LSTICC ), 1 ) 253: * 254: * NOTE: The computation of VN1( LSTICC ) relies on the fact that 255: * SNRM2 does not fail on vectors with norm below the value of 256: * SQRT(DLAMCH('S')) 257: * 258: VN2( LSTICC ) = VN1( LSTICC ) 259: LSTICC = ITEMP 260: GO TO 60 261: END IF 262: * 263: RETURN 264: * 265: * End of ZLAQPS 266: * 267: END