![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) 2: * 3: * -- LAPACK auxiliary routine (version 3.2.1) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * -- April 2009 -- 7: * 8: * .. Scalar Arguments .. 9: INTEGER K, LDA, LDT, LDY, N, NB 10: * .. 11: * .. Array Arguments .. 12: COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ), 13: $ Y( LDY, NB ) 14: * .. 15: * 16: * Purpose 17: * ======= 18: * 19: * ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1) 20: * matrix A so that elements below the k-th subdiagonal are zero. The 21: * reduction is performed by an unitary similarity transformation 22: * Q' * A * Q. The routine returns the matrices V and T which determine 23: * Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. 24: * 25: * This is an auxiliary routine called by ZGEHRD. 26: * 27: * Arguments 28: * ========= 29: * 30: * N (input) INTEGER 31: * The order of the matrix A. 32: * 33: * K (input) INTEGER 34: * The offset for the reduction. Elements below the k-th 35: * subdiagonal in the first NB columns are reduced to zero. 36: * K < N. 37: * 38: * NB (input) INTEGER 39: * The number of columns to be reduced. 40: * 41: * A (input/output) COMPLEX*16 array, dimension (LDA,N-K+1) 42: * On entry, the n-by-(n-k+1) general matrix A. 43: * On exit, the elements on and above the k-th subdiagonal in 44: * the first NB columns are overwritten with the corresponding 45: * elements of the reduced matrix; the elements below the k-th 46: * subdiagonal, with the array TAU, represent the matrix Q as a 47: * product of elementary reflectors. The other columns of A are 48: * unchanged. See Further Details. 49: * 50: * LDA (input) INTEGER 51: * The leading dimension of the array A. LDA >= max(1,N). 52: * 53: * TAU (output) COMPLEX*16 array, dimension (NB) 54: * The scalar factors of the elementary reflectors. See Further 55: * Details. 56: * 57: * T (output) COMPLEX*16 array, dimension (LDT,NB) 58: * The upper triangular matrix T. 59: * 60: * LDT (input) INTEGER 61: * The leading dimension of the array T. LDT >= NB. 62: * 63: * Y (output) COMPLEX*16 array, dimension (LDY,NB) 64: * The n-by-nb matrix Y. 65: * 66: * LDY (input) INTEGER 67: * The leading dimension of the array Y. LDY >= N. 68: * 69: * Further Details 70: * =============== 71: * 72: * The matrix Q is represented as a product of nb elementary reflectors 73: * 74: * Q = H(1) H(2) . . . H(nb). 75: * 76: * Each H(i) has the form 77: * 78: * H(i) = I - tau * v * v' 79: * 80: * where tau is a complex scalar, and v is a complex vector with 81: * v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in 82: * A(i+k+1:n,i), and tau in TAU(i). 83: * 84: * The elements of the vectors v together form the (n-k+1)-by-nb matrix 85: * V which is needed, with T and Y, to apply the transformation to the 86: * unreduced part of the matrix, using an update of the form: 87: * A := (I - V*T*V') * (A - Y*V'). 88: * 89: * The contents of A on exit are illustrated by the following example 90: * with n = 7, k = 3 and nb = 2: 91: * 92: * ( a a a a a ) 93: * ( a a a a a ) 94: * ( a a a a a ) 95: * ( h h a a a ) 96: * ( v1 h a a a ) 97: * ( v1 v2 a a a ) 98: * ( v1 v2 a a a ) 99: * 100: * where a denotes an element of the original matrix A, h denotes a 101: * modified element of the upper Hessenberg matrix H, and vi denotes an 102: * element of the vector defining H(i). 103: * 104: * This subroutine is a slight modification of LAPACK-3.0's DLAHRD 105: * incorporating improvements proposed by Quintana-Orti and Van de 106: * Gejin. Note that the entries of A(1:K,2:NB) differ from those 107: * returned by the original LAPACK-3.0's DLAHRD routine. (This 108: * subroutine is not backward compatible with LAPACK-3.0's DLAHRD.) 109: * 110: * References 111: * ========== 112: * 113: * Gregorio Quintana-Orti and Robert van de Geijn, "Improving the 114: * performance of reduction to Hessenberg form," ACM Transactions on 115: * Mathematical Software, 32(2):180-194, June 2006. 116: * 117: * ===================================================================== 118: * 119: * .. Parameters .. 120: COMPLEX*16 ZERO, ONE 121: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 122: $ ONE = ( 1.0D+0, 0.0D+0 ) ) 123: * .. 124: * .. Local Scalars .. 125: INTEGER I 126: COMPLEX*16 EI 127: * .. 128: * .. External Subroutines .. 129: EXTERNAL ZAXPY, ZCOPY, ZGEMM, ZGEMV, ZLACPY, 130: $ ZLARFG, ZSCAL, ZTRMM, ZTRMV, ZLACGV 131: * .. 132: * .. Intrinsic Functions .. 133: INTRINSIC MIN 134: * .. 135: * .. Executable Statements .. 136: * 137: * Quick return if possible 138: * 139: IF( N.LE.1 ) 140: $ RETURN 141: * 142: DO 10 I = 1, NB 143: IF( I.GT.1 ) THEN 144: * 145: * Update A(K+1:N,I) 146: * 147: * Update I-th column of A - Y * V' 148: * 149: CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) 150: CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, 151: $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) 152: CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) 153: * 154: * Apply I - V * T' * V' to this column (call it b) from the 155: * left, using the last column of T as workspace 156: * 157: * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) 158: * ( V2 ) ( b2 ) 159: * 160: * where V1 is unit lower triangular 161: * 162: * w := V1' * b1 163: * 164: CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) 165: CALL ZTRMV( 'Lower', 'Conjugate transpose', 'UNIT', 166: $ I-1, A( K+1, 1 ), 167: $ LDA, T( 1, NB ), 1 ) 168: * 169: * w := w + V2'*b2 170: * 171: CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, 172: $ ONE, A( K+I, 1 ), 173: $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) 174: * 175: * w := T'*w 176: * 177: CALL ZTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT', 178: $ I-1, T, LDT, 179: $ T( 1, NB ), 1 ) 180: * 181: * b2 := b2 - V2*w 182: * 183: CALL ZGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, 184: $ A( K+I, 1 ), 185: $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) 186: * 187: * b1 := b1 - V1*w 188: * 189: CALL ZTRMV( 'Lower', 'NO TRANSPOSE', 190: $ 'UNIT', I-1, 191: $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) 192: CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) 193: * 194: A( K+I-1, I-1 ) = EI 195: END IF 196: * 197: * Generate the elementary reflector H(I) to annihilate 198: * A(K+I+1:N,I) 199: * 200: CALL ZLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, 201: $ TAU( I ) ) 202: EI = A( K+I, I ) 203: A( K+I, I ) = ONE 204: * 205: * Compute Y(K+1:N,I) 206: * 207: CALL ZGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, 208: $ ONE, A( K+1, I+1 ), 209: $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) 210: CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, 211: $ ONE, A( K+I, 1 ), LDA, 212: $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) 213: CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, 214: $ Y( K+1, 1 ), LDY, 215: $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) 216: CALL ZSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) 217: * 218: * Compute T(1:I,I) 219: * 220: CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) 221: CALL ZTRMV( 'Upper', 'No Transpose', 'NON-UNIT', 222: $ I-1, T, LDT, 223: $ T( 1, I ), 1 ) 224: T( I, I ) = TAU( I ) 225: * 226: 10 CONTINUE 227: A( K+NB, NB ) = EI 228: * 229: * Compute Y(1:K,1:NB) 230: * 231: CALL ZLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) 232: CALL ZTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', 233: $ 'UNIT', K, NB, 234: $ ONE, A( K+1, 1 ), LDA, Y, LDY ) 235: IF( N.GT.K+NB ) 236: $ CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, 237: $ NB, N-K-NB, ONE, 238: $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, 239: $ LDY ) 240: CALL ZTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', 241: $ 'NON-UNIT', K, NB, 242: $ ONE, T, LDT, Y, LDY ) 243: * 244: RETURN 245: * 246: * End of ZLAHR2 247: * 248: END