![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
1: SUBROUTINE DLAQPS( 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 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ), 15: $ VN1( * ), VN2( * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * DLAQPS computes a step of QR factorization with column pivoting 22: * of a real 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (NB) 76: * Auxiliar vector. 77: * 78: * F (input/output) DOUBLE PRECISION 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: * Partial column norm updating strategy modified by 92: * Z. Drmac and Z. Bujanovic, Dept. of Mathematics, 93: * University of Zagreb, Croatia. 94: * June 2010 95: * For more details see LAPACK Working Note 176. 96: * ===================================================================== 97: * 98: * .. Parameters .. 99: DOUBLE PRECISION ZERO, ONE 100: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 101: * .. 102: * .. Local Scalars .. 103: INTEGER ITEMP, J, K, LASTRK, LSTICC, PVT, RK 104: DOUBLE PRECISION AKK, TEMP, TEMP2, TOL3Z 105: * .. 106: * .. External Subroutines .. 107: EXTERNAL DGEMM, DGEMV, DLARFG, DSWAP 108: * .. 109: * .. Intrinsic Functions .. 110: INTRINSIC ABS, DBLE, MAX, MIN, NINT, SQRT 111: * .. 112: * .. External Functions .. 113: INTEGER IDAMAX 114: DOUBLE PRECISION DLAMCH, DNRM2 115: EXTERNAL IDAMAX, DLAMCH, DNRM2 116: * .. 117: * .. Executable Statements .. 118: * 119: LASTRK = MIN( M, N+OFFSET ) 120: LSTICC = 0 121: K = 0 122: TOL3Z = SQRT(DLAMCH('Epsilon')) 123: * 124: * Beginning of while loop. 125: * 126: 10 CONTINUE 127: IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN 128: K = K + 1 129: RK = OFFSET + K 130: * 131: * Determine ith pivot column and swap if necessary 132: * 133: PVT = ( K-1 ) + IDAMAX( N-K+1, VN1( K ), 1 ) 134: IF( PVT.NE.K ) THEN 135: CALL DSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 ) 136: CALL DSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF ) 137: ITEMP = JPVT( PVT ) 138: JPVT( PVT ) = JPVT( K ) 139: JPVT( K ) = ITEMP 140: VN1( PVT ) = VN1( K ) 141: VN2( PVT ) = VN2( K ) 142: END IF 143: * 144: * Apply previous Householder reflectors to column K: 145: * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'. 146: * 147: IF( K.GT.1 ) THEN 148: CALL DGEMV( 'No transpose', M-RK+1, K-1, -ONE, A( RK, 1 ), 149: $ LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 ) 150: END IF 151: * 152: * Generate elementary reflector H(k). 153: * 154: IF( RK.LT.M ) THEN 155: CALL DLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) 156: ELSE 157: CALL DLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) 158: END IF 159: * 160: AKK = A( RK, K ) 161: A( RK, K ) = ONE 162: * 163: * Compute Kth column of F: 164: * 165: * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K). 166: * 167: IF( K.LT.N ) THEN 168: CALL DGEMV( 'Transpose', M-RK+1, N-K, TAU( K ), 169: $ A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO, 170: $ F( K+1, K ), 1 ) 171: END IF 172: * 173: * Padding F(1:K,K) with zeros. 174: * 175: DO 20 J = 1, K 176: F( J, K ) = ZERO 177: 20 CONTINUE 178: * 179: * Incremental updating of F: 180: * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)' 181: * *A(RK:M,K). 182: * 183: IF( K.GT.1 ) THEN 184: CALL DGEMV( 'Transpose', M-RK+1, K-1, -TAU( K ), A( RK, 1 ), 185: $ LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 ) 186: * 187: CALL DGEMV( 'No transpose', N, K-1, ONE, F( 1, 1 ), LDF, 188: $ AUXV( 1 ), 1, ONE, F( 1, K ), 1 ) 189: END IF 190: * 191: * Update the current row of A: 192: * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'. 193: * 194: IF( K.LT.N ) THEN 195: CALL DGEMV( 'No transpose', N-K, K, -ONE, F( K+1, 1 ), LDF, 196: $ A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA ) 197: END IF 198: * 199: * Update partial column norms. 200: * 201: IF( RK.LT.LASTRK ) THEN 202: DO 30 J = K + 1, N 203: IF( VN1( J ).NE.ZERO ) THEN 204: * 205: * NOTE: The following 4 lines follow from the analysis in 206: * Lapack Working Note 176. 207: * 208: TEMP = ABS( A( RK, J ) ) / VN1( J ) 209: TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) 210: TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 211: IF( TEMP2 .LE. TOL3Z ) THEN 212: VN2( J ) = DBLE( LSTICC ) 213: LSTICC = J 214: ELSE 215: VN1( J ) = VN1( J )*SQRT( TEMP ) 216: END IF 217: END IF 218: 30 CONTINUE 219: END IF 220: * 221: A( RK, K ) = AKK 222: * 223: * End of while loop. 224: * 225: GO TO 10 226: END IF 227: KB = K 228: RK = OFFSET + KB 229: * 230: * Apply the block reflector to the rest of the matrix: 231: * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - 232: * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'. 233: * 234: IF( KB.LT.MIN( N, M-OFFSET ) ) THEN 235: CALL DGEMM( 'No transpose', 'Transpose', M-RK, N-KB, KB, -ONE, 236: $ A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE, 237: $ A( RK+1, KB+1 ), LDA ) 238: END IF 239: * 240: * Recomputation of difficult columns. 241: * 242: 40 CONTINUE 243: IF( LSTICC.GT.0 ) THEN 244: ITEMP = NINT( VN2( LSTICC ) ) 245: VN1( LSTICC ) = DNRM2( M-RK, A( RK+1, LSTICC ), 1 ) 246: * 247: * NOTE: The computation of VN1( LSTICC ) relies on the fact that 248: * SNRM2 does not fail on vectors with norm below the value of 249: * SQRT(DLAMCH('S')) 250: * 251: VN2( LSTICC ) = VN1( LSTICC ) 252: LSTICC = ITEMP 253: GO TO 40 254: END IF 255: * 256: RETURN 257: * 258: * End of DLAQPS 259: * 260: END