Annotation of rpl/lapack/lapack/dgehrd.f, revision 1.7
1.1 bertrand 1: SUBROUTINE DGEHRD( 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: DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
13: * ..
14: *
15: * Purpose
16: * =======
17: *
18: * DGEHRD reduces a real general matrix A to upper Hessenberg form H by
19: * an orthogonal 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 DGEBAL; 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) DOUBLE PRECISION 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 orthogonal 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 real scalar, and v is a real 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: DOUBLE PRECISION ZERO, ONE
111: PARAMETER ( ZERO = 0.0D+0,
112: $ ONE = 1.0D+0 )
113: * ..
114: * .. Local Scalars ..
115: LOGICAL LQUERY
116: INTEGER I, IB, IINFO, IWS, J, LDWORK, LWKOPT, NB,
117: $ NBMIN, NH, NX
118: DOUBLE PRECISION EI
119: * ..
120: * .. Local Arrays ..
121: DOUBLE PRECISION T( LDT, NBMAX )
122: * ..
123: * .. External Subroutines ..
124: EXTERNAL DAXPY, DGEHD2, DGEMM, DLAHR2, DLARFB, DTRMM,
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, 'DGEHRD', ' ', 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( 'DGEHRD', -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, 'DGEHRD', ' ', 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, 'DGEHRD', ' ', 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, 'DGEHRD', ' ', 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 DLAHR2( 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 DGEMM( 'No transpose', '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 DTRMM( 'Right', 'Lower', 'Transpose',
248: $ 'Unit', I, IB-1,
249: $ ONE, A( I+1, I ), LDA, WORK, LDWORK )
250: DO 30 J = 0, IB-2
251: CALL DAXPY( 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 DLARFB( 'Left', '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 DGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
268: WORK( 1 ) = IWS
269: *
270: RETURN
271: *
272: * End of DGEHRD
273: *
274: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgehrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack