![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 2: IMPLICIT NONE 3: * 4: * -- LAPACK auxiliary routine (version 3.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: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER DIRECT, STOREV 11: INTEGER K, LDT, LDV, N 12: * .. 13: * .. Array Arguments .. 14: DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLARFT forms the triangular factor T of a real block reflector H 21: * of order n, which is defined as a product of k elementary reflectors. 22: * 23: * If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; 24: * 25: * If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. 26: * 27: * If STOREV = 'C', the vector which defines the elementary reflector 28: * H(i) is stored in the i-th column of the array V, and 29: * 30: * H = I - V * T * V' 31: * 32: * If STOREV = 'R', the vector which defines the elementary reflector 33: * H(i) is stored in the i-th row of the array V, and 34: * 35: * H = I - V' * T * V 36: * 37: * Arguments 38: * ========= 39: * 40: * DIRECT (input) CHARACTER*1 41: * Specifies the order in which the elementary reflectors are 42: * multiplied to form the block reflector: 43: * = 'F': H = H(1) H(2) . . . H(k) (Forward) 44: * = 'B': H = H(k) . . . H(2) H(1) (Backward) 45: * 46: * STOREV (input) CHARACTER*1 47: * Specifies how the vectors which define the elementary 48: * reflectors are stored (see also Further Details): 49: * = 'C': columnwise 50: * = 'R': rowwise 51: * 52: * N (input) INTEGER 53: * The order of the block reflector H. N >= 0. 54: * 55: * K (input) INTEGER 56: * The order of the triangular factor T (= the number of 57: * elementary reflectors). K >= 1. 58: * 59: * V (input/output) DOUBLE PRECISION array, dimension 60: * (LDV,K) if STOREV = 'C' 61: * (LDV,N) if STOREV = 'R' 62: * The matrix V. See further details. 63: * 64: * LDV (input) INTEGER 65: * The leading dimension of the array V. 66: * If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. 67: * 68: * TAU (input) DOUBLE PRECISION array, dimension (K) 69: * TAU(i) must contain the scalar factor of the elementary 70: * reflector H(i). 71: * 72: * T (output) DOUBLE PRECISION array, dimension (LDT,K) 73: * The k by k triangular factor T of the block reflector. 74: * If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is 75: * lower triangular. The rest of the array is not used. 76: * 77: * LDT (input) INTEGER 78: * The leading dimension of the array T. LDT >= K. 79: * 80: * Further Details 81: * =============== 82: * 83: * The shape of the matrix V and the storage of the vectors which define 84: * the H(i) is best illustrated by the following example with n = 5 and 85: * k = 3. The elements equal to 1 are not stored; the corresponding 86: * array elements are modified but restored on exit. The rest of the 87: * array is not used. 88: * 89: * DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 90: * 91: * V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 92: * ( v1 1 ) ( 1 v2 v2 v2 ) 93: * ( v1 v2 1 ) ( 1 v3 v3 ) 94: * ( v1 v2 v3 ) 95: * ( v1 v2 v3 ) 96: * 97: * DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 98: * 99: * V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 100: * ( v1 v2 v3 ) ( v2 v2 v2 1 ) 101: * ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 102: * ( 1 v3 ) 103: * ( 1 ) 104: * 105: * ===================================================================== 106: * 107: * .. Parameters .. 108: DOUBLE PRECISION ONE, ZERO 109: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 110: * .. 111: * .. Local Scalars .. 112: INTEGER I, J, PREVLASTV, LASTV 113: DOUBLE PRECISION VII 114: * .. 115: * .. External Subroutines .. 116: EXTERNAL DGEMV, DTRMV 117: * .. 118: * .. External Functions .. 119: LOGICAL LSAME 120: EXTERNAL LSAME 121: * .. 122: * .. Executable Statements .. 123: * 124: * Quick return if possible 125: * 126: IF( N.EQ.0 ) 127: $ RETURN 128: * 129: IF( LSAME( DIRECT, 'F' ) ) THEN 130: PREVLASTV = N 131: DO 20 I = 1, K 132: PREVLASTV = MAX( I, PREVLASTV ) 133: IF( TAU( I ).EQ.ZERO ) THEN 134: * 135: * H(i) = I 136: * 137: DO 10 J = 1, I 138: T( J, I ) = ZERO 139: 10 CONTINUE 140: ELSE 141: * 142: * general case 143: * 144: VII = V( I, I ) 145: V( I, I ) = ONE 146: IF( LSAME( STOREV, 'C' ) ) THEN 147: ! Skip any trailing zeros. 148: DO LASTV = N, I+1, -1 149: IF( V( LASTV, I ).NE.ZERO ) EXIT 150: END DO 151: J = MIN( LASTV, PREVLASTV ) 152: * 153: * T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)' * V(i:j,i) 154: * 155: CALL DGEMV( 'Transpose', J-I+1, I-1, -TAU( I ), 156: $ V( I, 1 ), LDV, V( I, I ), 1, ZERO, 157: $ T( 1, I ), 1 ) 158: ELSE 159: ! Skip any trailing zeros. 160: DO LASTV = N, I+1, -1 161: IF( V( I, LASTV ).NE.ZERO ) EXIT 162: END DO 163: J = MIN( LASTV, PREVLASTV ) 164: * 165: * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)' 166: * 167: CALL DGEMV( 'No transpose', I-1, J-I+1, -TAU( I ), 168: $ V( 1, I ), LDV, V( I, I ), LDV, ZERO, 169: $ T( 1, I ), 1 ) 170: END IF 171: V( I, I ) = VII 172: * 173: * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 174: * 175: CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 176: $ LDT, T( 1, I ), 1 ) 177: T( I, I ) = TAU( I ) 178: IF( I.GT.1 ) THEN 179: PREVLASTV = MAX( PREVLASTV, LASTV ) 180: ELSE 181: PREVLASTV = LASTV 182: END IF 183: END IF 184: 20 CONTINUE 185: ELSE 186: PREVLASTV = 1 187: DO 40 I = K, 1, -1 188: IF( TAU( I ).EQ.ZERO ) THEN 189: * 190: * H(i) = I 191: * 192: DO 30 J = I, K 193: T( J, I ) = ZERO 194: 30 CONTINUE 195: ELSE 196: * 197: * general case 198: * 199: IF( I.LT.K ) THEN 200: IF( LSAME( STOREV, 'C' ) ) THEN 201: VII = V( N-K+I, I ) 202: V( N-K+I, I ) = ONE 203: ! Skip any leading zeros. 204: DO LASTV = 1, I-1 205: IF( V( LASTV, I ).NE.ZERO ) EXIT 206: END DO 207: J = MAX( LASTV, PREVLASTV ) 208: * 209: * T(i+1:k,i) := 210: * - tau(i) * V(j:n-k+i,i+1:k)' * V(j:n-k+i,i) 211: * 212: CALL DGEMV( 'Transpose', N-K+I-J+1, K-I, -TAU( I ), 213: $ V( J, I+1 ), LDV, V( J, I ), 1, ZERO, 214: $ T( I+1, I ), 1 ) 215: V( N-K+I, I ) = VII 216: ELSE 217: VII = V( I, N-K+I ) 218: V( I, N-K+I ) = ONE 219: ! Skip any leading zeros. 220: DO LASTV = 1, I-1 221: IF( V( I, LASTV ).NE.ZERO ) EXIT 222: END DO 223: J = MAX( LASTV, PREVLASTV ) 224: * 225: * T(i+1:k,i) := 226: * - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)' 227: * 228: CALL DGEMV( 'No transpose', K-I, N-K+I-J+1, 229: $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, 230: $ ZERO, T( I+1, I ), 1 ) 231: V( I, N-K+I ) = VII 232: END IF 233: * 234: * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 235: * 236: CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 237: $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 238: IF( I.GT.1 ) THEN 239: PREVLASTV = MIN( PREVLASTV, LASTV ) 240: ELSE 241: PREVLASTV = LASTV 242: END IF 243: END IF 244: T( I, I ) = TAU( I ) 245: END IF 246: 40 CONTINUE 247: END IF 248: RETURN 249: * 250: * End of DLARFT 251: * 252: END