![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) 2: * 3: * -- LAPACK 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 IHI, ILO, INFO, LDA, LWORK, N 10: * .. 11: * .. Array Arguments .. 12: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) 13: * .. 14: * 15: * Purpose 16: * ======= 17: * 18: * ZGEHRD reduces a complex general matrix A to upper Hessenberg form H by 19: * an unitary similarity transformation: Q' * A * Q = H . 20: * 21: * Arguments 22: * ========= 23: * 24: * N (input) INTEGER 25: * The order of the matrix A. N >= 0. 26: * 27: * ILO (input) INTEGER 28: * IHI (input) INTEGER 29: * It is assumed that A is already upper triangular in rows 30: * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally 31: * set by a previous call to ZGEBAL; otherwise they should be 32: * set to 1 and N respectively. See Further Details. 33: * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 34: * 35: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 36: * On entry, the N-by-N general matrix to be reduced. 37: * On exit, the upper triangle and the first subdiagonal of A 38: * are overwritten with the upper Hessenberg matrix H, and the 39: * elements below the first subdiagonal, with the array TAU, 40: * represent the unitary matrix Q as a product of elementary 41: * reflectors. See Further Details. 42: * 43: * LDA (input) INTEGER 44: * The leading dimension of the array A. LDA >= max(1,N). 45: * 46: * TAU (output) COMPLEX*16 array, dimension (N-1) 47: * The scalar factors of the elementary reflectors (see Further 48: * Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to 49: * zero. 50: * 51: * WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) 52: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 53: * 54: * LWORK (input) INTEGER 55: * The length of the array WORK. LWORK >= max(1,N). 56: * For optimum performance LWORK >= N*NB, where NB is the 57: * optimal blocksize. 58: * 59: * If LWORK = -1, then a workspace query is assumed; the routine 60: * only calculates the optimal size of the WORK array, returns 61: * this value as the first entry of the WORK array, and no error 62: * message related to LWORK is issued by XERBLA. 63: * 64: * INFO (output) INTEGER 65: * = 0: successful exit 66: * < 0: if INFO = -i, the i-th argument had an illegal value. 67: * 68: * Further Details 69: * =============== 70: * 71: * The matrix Q is represented as a product of (ihi-ilo) elementary 72: * reflectors 73: * 74: * Q = H(ilo) H(ilo+1) . . . H(ihi-1). 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) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on 82: * exit in A(i+2:ihi,i), and tau in TAU(i). 83: * 84: * The contents of A are illustrated by the following example, with 85: * n = 7, ilo = 2 and ihi = 6: 86: * 87: * on entry, on exit, 88: * 89: * ( a a a a a a a ) ( a a h h h h a ) 90: * ( a a a a a a ) ( a h h h h a ) 91: * ( a a a a a a ) ( h h h h h h ) 92: * ( a a a a a a ) ( v2 h h h h h ) 93: * ( a a a a a a ) ( v2 v3 h h h h ) 94: * ( a a a a a a ) ( v2 v3 v4 h h h ) 95: * ( a ) ( a ) 96: * 97: * where a denotes an element of the original matrix A, h denotes a 98: * modified element of the upper Hessenberg matrix H, and vi denotes an 99: * element of the vector defining H(i). 100: * 101: * This file is a slight modification of LAPACK-3.0's DGEHRD 102: * subroutine incorporating improvements proposed by Quintana-Orti and 103: * Van de Geijn (2006). (See DLAHR2.) 104: * 105: * ===================================================================== 106: * 107: * .. Parameters .. 108: INTEGER NBMAX, LDT 109: PARAMETER ( NBMAX = 64, LDT = NBMAX+1 ) 110: COMPLEX*16 ZERO, ONE 111: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 112: $ ONE = ( 1.0D+0, 0.0D+0 ) ) 113: * .. 114: * .. Local Scalars .. 115: LOGICAL LQUERY 116: INTEGER I, IB, IINFO, IWS, J, LDWORK, LWKOPT, NB, 117: $ NBMIN, NH, NX 118: COMPLEX*16 EI 119: * .. 120: * .. Local Arrays .. 121: COMPLEX*16 T( LDT, NBMAX ) 122: * .. 123: * .. External Subroutines .. 124: EXTERNAL ZAXPY, ZGEHD2, ZGEMM, ZLAHR2, ZLARFB, ZTRMM, 125: $ XERBLA 126: * .. 127: * .. Intrinsic Functions .. 128: INTRINSIC MAX, MIN 129: * .. 130: * .. External Functions .. 131: INTEGER ILAENV 132: EXTERNAL ILAENV 133: * .. 134: * .. Executable Statements .. 135: * 136: * Test the input parameters 137: * 138: INFO = 0 139: NB = MIN( NBMAX, ILAENV( 1, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) ) 140: LWKOPT = N*NB 141: WORK( 1 ) = LWKOPT 142: LQUERY = ( LWORK.EQ.-1 ) 143: IF( N.LT.0 ) THEN 144: INFO = -1 145: ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN 146: INFO = -2 147: ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN 148: INFO = -3 149: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 150: INFO = -5 151: ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 152: INFO = -8 153: END IF 154: IF( INFO.NE.0 ) THEN 155: CALL XERBLA( 'ZGEHRD', -INFO ) 156: RETURN 157: ELSE IF( LQUERY ) THEN 158: RETURN 159: END IF 160: * 161: * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero 162: * 163: DO 10 I = 1, ILO - 1 164: TAU( I ) = ZERO 165: 10 CONTINUE 166: DO 20 I = MAX( 1, IHI ), N - 1 167: TAU( I ) = ZERO 168: 20 CONTINUE 169: * 170: * Quick return if possible 171: * 172: NH = IHI - ILO + 1 173: IF( NH.LE.1 ) THEN 174: WORK( 1 ) = 1 175: RETURN 176: END IF 177: * 178: * Determine the block size 179: * 180: NB = MIN( NBMAX, ILAENV( 1, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) ) 181: NBMIN = 2 182: IWS = 1 183: IF( NB.GT.1 .AND. NB.LT.NH ) THEN 184: * 185: * Determine when to cross over from blocked to unblocked code 186: * (last block is always handled by unblocked code) 187: * 188: NX = MAX( NB, ILAENV( 3, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) ) 189: IF( NX.LT.NH ) THEN 190: * 191: * Determine if workspace is large enough for blocked code 192: * 193: IWS = N*NB 194: IF( LWORK.LT.IWS ) THEN 195: * 196: * Not enough workspace to use optimal NB: determine the 197: * minimum value of NB, and reduce NB or force use of 198: * unblocked code 199: * 200: NBMIN = MAX( 2, ILAENV( 2, 'ZGEHRD', ' ', N, ILO, IHI, 201: $ -1 ) ) 202: IF( LWORK.GE.N*NBMIN ) THEN 203: NB = LWORK / N 204: ELSE 205: NB = 1 206: END IF 207: END IF 208: END IF 209: END IF 210: LDWORK = N 211: * 212: IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN 213: * 214: * Use unblocked code below 215: * 216: I = ILO 217: * 218: ELSE 219: * 220: * Use blocked code 221: * 222: DO 40 I = ILO, IHI - 1 - NX, NB 223: IB = MIN( NB, IHI-I ) 224: * 225: * Reduce columns i:i+ib-1 to Hessenberg form, returning the 226: * matrices V and T of the block reflector H = I - V*T*V' 227: * which performs the reduction, and also the matrix Y = A*V*T 228: * 229: CALL ZLAHR2( IHI, I, IB, A( 1, I ), LDA, TAU( I ), T, LDT, 230: $ WORK, LDWORK ) 231: * 232: * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the 233: * right, computing A := A - Y * V'. V(i+ib,ib-1) must be set 234: * to 1 235: * 236: EI = A( I+IB, I+IB-1 ) 237: A( I+IB, I+IB-1 ) = ONE 238: CALL ZGEMM( 'No transpose', 'Conjugate transpose', 239: $ IHI, IHI-I-IB+1, 240: $ IB, -ONE, WORK, LDWORK, A( I+IB, I ), LDA, ONE, 241: $ A( 1, I+IB ), LDA ) 242: A( I+IB, I+IB-1 ) = EI 243: * 244: * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the 245: * right 246: * 247: CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose', 248: $ 'Unit', I, IB-1, 249: $ ONE, A( I+1, I ), LDA, WORK, LDWORK ) 250: DO 30 J = 0, IB-2 251: CALL ZAXPY( I, -ONE, WORK( LDWORK*J+1 ), 1, 252: $ A( 1, I+J+1 ), 1 ) 253: 30 CONTINUE 254: * 255: * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the 256: * left 257: * 258: CALL ZLARFB( 'Left', 'Conjugate transpose', 'Forward', 259: $ 'Columnwise', 260: $ IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA, T, LDT, 261: $ A( I+1, I+IB ), LDA, WORK, LDWORK ) 262: 40 CONTINUE 263: END IF 264: * 265: * Use unblocked code to reduce the rest of the matrix 266: * 267: CALL ZGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO ) 268: WORK( 1 ) = IWS 269: * 270: RETURN 271: * 272: * End of ZGEHRD 273: * 274: END