![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO ) 2: * 3: * -- LAPACK deprecated computational routine (version 3.2.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * June 2010 7: * 8: * .. Scalar Arguments .. 9: INTEGER INFO, LDA, M, N 10: * .. 11: * .. Array Arguments .. 12: INTEGER JPVT( * ) 13: DOUBLE PRECISION RWORK( * ) 14: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * This routine is deprecated and has been replaced by routine ZGEQP3. 21: * 22: * ZGEQPF computes a QR factorization with column pivoting of a 23: * complex M-by-N matrix A: A*P = Q*R. 24: * 25: * Arguments 26: * ========= 27: * 28: * M (input) INTEGER 29: * The number of rows of the matrix A. M >= 0. 30: * 31: * N (input) INTEGER 32: * The number of columns of the matrix A. N >= 0 33: * 34: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 35: * On entry, the M-by-N matrix A. 36: * On exit, the upper triangle of the array contains the 37: * min(M,N)-by-N upper triangular matrix R; the elements 38: * below the diagonal, together with the array TAU, 39: * represent the unitary matrix Q as a product of 40: * min(m,n) elementary reflectors. 41: * 42: * LDA (input) INTEGER 43: * The leading dimension of the array A. LDA >= max(1,M). 44: * 45: * JPVT (input/output) INTEGER array, dimension (N) 46: * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted 47: * to the front of A*P (a leading column); if JPVT(i) = 0, 48: * the i-th column of A is a free column. 49: * On exit, if JPVT(i) = k, then the i-th column of A*P 50: * was the k-th column of A. 51: * 52: * TAU (output) COMPLEX*16 array, dimension (min(M,N)) 53: * The scalar factors of the elementary reflectors. 54: * 55: * WORK (workspace) COMPLEX*16 array, dimension (N) 56: * 57: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 58: * 59: * INFO (output) INTEGER 60: * = 0: successful exit 61: * < 0: if INFO = -i, the i-th argument had an illegal value 62: * 63: * Further Details 64: * =============== 65: * 66: * The matrix Q is represented as a product of elementary reflectors 67: * 68: * Q = H(1) H(2) . . . H(n) 69: * 70: * Each H(i) has the form 71: * 72: * H = I - tau * v * v' 73: * 74: * where tau is a complex scalar, and v is a complex vector with 75: * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). 76: * 77: * The matrix P is represented in jpvt as follows: If 78: * jpvt(j) = i 79: * then the jth column of P is the ith canonical unit vector. 80: * 81: * Partial column norm updating strategy modified by 82: * Z. Drmac and Z. Bujanovic, Dept. of Mathematics, 83: * University of Zagreb, Croatia. 84: * June 2010 85: * For more details see LAPACK Working Note 176. 86: * 87: * ===================================================================== 88: * 89: * .. Parameters .. 90: DOUBLE PRECISION ZERO, ONE 91: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 92: * .. 93: * .. Local Scalars .. 94: INTEGER I, ITEMP, J, MA, MN, PVT 95: DOUBLE PRECISION TEMP, TEMP2, TOL3Z 96: COMPLEX*16 AII 97: * .. 98: * .. External Subroutines .. 99: EXTERNAL XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R 100: * .. 101: * .. Intrinsic Functions .. 102: INTRINSIC ABS, DCMPLX, DCONJG, MAX, MIN, SQRT 103: * .. 104: * .. External Functions .. 105: INTEGER IDAMAX 106: DOUBLE PRECISION DLAMCH, DZNRM2 107: EXTERNAL IDAMAX, DLAMCH, DZNRM2 108: * .. 109: * .. Executable Statements .. 110: * 111: * Test the input arguments 112: * 113: INFO = 0 114: IF( M.LT.0 ) THEN 115: INFO = -1 116: ELSE IF( N.LT.0 ) THEN 117: INFO = -2 118: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 119: INFO = -4 120: END IF 121: IF( INFO.NE.0 ) THEN 122: CALL XERBLA( 'ZGEQPF', -INFO ) 123: RETURN 124: END IF 125: * 126: MN = MIN( M, N ) 127: TOL3Z = SQRT(DLAMCH('Epsilon')) 128: * 129: * Move initial columns up front 130: * 131: ITEMP = 1 132: DO 10 I = 1, N 133: IF( JPVT( I ).NE.0 ) THEN 134: IF( I.NE.ITEMP ) THEN 135: CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 ) 136: JPVT( I ) = JPVT( ITEMP ) 137: JPVT( ITEMP ) = I 138: ELSE 139: JPVT( I ) = I 140: END IF 141: ITEMP = ITEMP + 1 142: ELSE 143: JPVT( I ) = I 144: END IF 145: 10 CONTINUE 146: ITEMP = ITEMP - 1 147: * 148: * Compute the QR factorization and update remaining columns 149: * 150: IF( ITEMP.GT.0 ) THEN 151: MA = MIN( ITEMP, M ) 152: CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO ) 153: IF( MA.LT.N ) THEN 154: CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A, 155: $ LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO ) 156: END IF 157: END IF 158: * 159: IF( ITEMP.LT.MN ) THEN 160: * 161: * Initialize partial column norms. The first n elements of 162: * work store the exact column norms. 163: * 164: DO 20 I = ITEMP + 1, N 165: RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 ) 166: RWORK( N+I ) = RWORK( I ) 167: 20 CONTINUE 168: * 169: * Compute factorization 170: * 171: DO 40 I = ITEMP + 1, MN 172: * 173: * Determine ith pivot column and swap if necessary 174: * 175: PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 ) 176: * 177: IF( PVT.NE.I ) THEN 178: CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) 179: ITEMP = JPVT( PVT ) 180: JPVT( PVT ) = JPVT( I ) 181: JPVT( I ) = ITEMP 182: RWORK( PVT ) = RWORK( I ) 183: RWORK( N+PVT ) = RWORK( N+I ) 184: END IF 185: * 186: * Generate elementary reflector H(i) 187: * 188: AII = A( I, I ) 189: CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1, 190: $ TAU( I ) ) 191: A( I, I ) = AII 192: * 193: IF( I.LT.N ) THEN 194: * 195: * Apply H(i) to A(i:m,i+1:n) from the left 196: * 197: AII = A( I, I ) 198: A( I, I ) = DCMPLX( ONE ) 199: CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, 200: $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) 201: A( I, I ) = AII 202: END IF 203: * 204: * Update partial column norms 205: * 206: DO 30 J = I + 1, N 207: IF( RWORK( J ).NE.ZERO ) THEN 208: * 209: * NOTE: The following 4 lines follow from the analysis in 210: * Lapack Working Note 176. 211: * 212: TEMP = ABS( A( I, J ) ) / RWORK( J ) 213: TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) ) 214: TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2 215: IF( TEMP2 .LE. TOL3Z ) THEN 216: IF( M-I.GT.0 ) THEN 217: RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 ) 218: RWORK( N+J ) = RWORK( J ) 219: ELSE 220: RWORK( J ) = ZERO 221: RWORK( N+J ) = ZERO 222: END IF 223: ELSE 224: RWORK( J ) = RWORK( J )*SQRT( TEMP ) 225: END IF 226: END IF 227: 30 CONTINUE 228: * 229: 40 CONTINUE 230: END IF 231: RETURN 232: * 233: * End of ZGEQPF 234: * 235: END