LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgehrd.f
Go to the documentation of this file.
1*> \brief \b DGEHRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGEHRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgehrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgehrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgehrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER IHI, ILO, INFO, LDA, LWORK, N
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> DGEHRD reduces a real general matrix A to upper Hessenberg form H by
37*> an orthogonal similarity transformation: Q**T * A * Q = H .
38*> \endverbatim
39*
40* Arguments:
41* ==========
42*
43*> \param[in] N
44*> \verbatim
45*> N is INTEGER
46*> The order of the matrix A. N >= 0.
47*> \endverbatim
48*>
49*> \param[in] ILO
50*> \verbatim
51*> ILO is INTEGER
52*> \endverbatim
53*>
54*> \param[in] IHI
55*> \verbatim
56*> IHI is INTEGER
57*>
58*> It is assumed that A is already upper triangular in rows
59*> and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
60*> set by a previous call to DGEBAL; otherwise they should be
61*> set to 1 and N respectively. See Further Details.
62*> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
63*> \endverbatim
64*>
65*> \param[in,out] A
66*> \verbatim
67*> A is DOUBLE PRECISION array, dimension (LDA,N)
68*> On entry, the N-by-N general matrix to be reduced.
69*> On exit, the upper triangle and the first subdiagonal of A
70*> are overwritten with the upper Hessenberg matrix H, and the
71*> elements below the first subdiagonal, with the array TAU,
72*> represent the orthogonal matrix Q as a product of elementary
73*> reflectors. See Further Details.
74*> \endverbatim
75*>
76*> \param[in] LDA
77*> \verbatim
78*> LDA is INTEGER
79*> The leading dimension of the array A. LDA >= max(1,N).
80*> \endverbatim
81*>
82*> \param[out] TAU
83*> \verbatim
84*> TAU is DOUBLE PRECISION array, dimension (N-1)
85*> The scalar factors of the elementary reflectors (see Further
86*> Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
87*> zero.
88*> \endverbatim
89*>
90*> \param[out] WORK
91*> \verbatim
92*> WORK is DOUBLE PRECISION array, dimension (LWORK)
93*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
94*> \endverbatim
95*>
96*> \param[in] LWORK
97*> \verbatim
98*> LWORK is INTEGER
99*> The length of the array WORK. LWORK >= max(1,N).
100*> For good performance, LWORK should generally be larger.
101*>
102*> If LWORK = -1, then a workspace query is assumed; the routine
103*> only calculates the optimal size of the WORK array, returns
104*> this value as the first entry of the WORK array, and no error
105*> message related to LWORK is issued by XERBLA.
106*> \endverbatim
107*>
108*> \param[out] INFO
109*> \verbatim
110*> INFO is INTEGER
111*> = 0: successful exit
112*> < 0: if INFO = -i, the i-th argument had an illegal value.
113*> \endverbatim
114*
115* Authors:
116* ========
117*
118*> \author Univ. of Tennessee
119*> \author Univ. of California Berkeley
120*> \author Univ. of Colorado Denver
121*> \author NAG Ltd.
122*
123*> \ingroup gehrd
124*
125*> \par Further Details:
126* =====================
127*>
128*> \verbatim
129*>
130*> The matrix Q is represented as a product of (ihi-ilo) elementary
131*> reflectors
132*>
133*> Q = H(ilo) H(ilo+1) . . . H(ihi-1).
134*>
135*> Each H(i) has the form
136*>
137*> H(i) = I - tau * v * v**T
138*>
139*> where tau is a real scalar, and v is a real vector with
140*> v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
141*> exit in A(i+2:ihi,i), and tau in TAU(i).
142*>
143*> The contents of A are illustrated by the following example, with
144*> n = 7, ilo = 2 and ihi = 6:
145*>
146*> on entry, on exit,
147*>
148*> ( a a a a a a a ) ( a a h h h h a )
149*> ( a a a a a a ) ( a h h h h a )
150*> ( a a a a a a ) ( h h h h h h )
151*> ( a a a a a a ) ( v2 h h h h h )
152*> ( a a a a a a ) ( v2 v3 h h h h )
153*> ( a a a a a a ) ( v2 v3 v4 h h h )
154*> ( a ) ( a )
155*>
156*> where a denotes an element of the original matrix A, h denotes a
157*> modified element of the upper Hessenberg matrix H, and vi denotes an
158*> element of the vector defining H(i).
159*>
160*> This file is a slight modification of LAPACK-3.0's DGEHRD
161*> subroutine incorporating improvements proposed by Quintana-Orti and
162*> Van de Geijn (2006). (See DLAHR2.)
163*> \endverbatim
164*>
165* =====================================================================
166 SUBROUTINE dgehrd( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
167*
168* -- LAPACK computational routine --
169* -- LAPACK is a software package provided by Univ. of Tennessee, --
170* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171*
172* .. Scalar Arguments ..
173 INTEGER IHI, ILO, INFO, LDA, LWORK, N
174* ..
175* .. Array Arguments ..
176 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 INTEGER NBMAX, LDT, TSIZE
183 parameter( nbmax = 64, ldt = nbmax+1,
184 $ tsize = ldt*nbmax )
185 DOUBLE PRECISION ZERO, ONE
186 parameter( zero = 0.0d+0,
187 $ one = 1.0d+0 )
188* ..
189* .. Local Scalars ..
190 LOGICAL LQUERY
191 INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
192 $ NBMIN, NH, NX
193 DOUBLE PRECISION EI
194* ..
195* .. External Subroutines ..
196 EXTERNAL daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm,
197 $ xerbla
198* ..
199* .. Intrinsic Functions ..
200 INTRINSIC max, min
201* ..
202* .. External Functions ..
203 INTEGER ILAENV
204 EXTERNAL ilaenv
205* ..
206* .. Executable Statements ..
207*
208* Test the input parameters
209*
210 info = 0
211 lquery = ( lwork.EQ.-1 )
212 IF( n.LT.0 ) THEN
213 info = -1
214 ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
215 info = -2
216 ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
217 info = -3
218 ELSE IF( lda.LT.max( 1, n ) ) THEN
219 info = -5
220 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
221 info = -8
222 END IF
223*
224 IF( info.EQ.0 ) THEN
225*
226* Compute the workspace requirements
227*
228 nb = min( nbmax, ilaenv( 1, 'DGEHRD', ' ', n, ilo, ihi, -1 ) )
229 lwkopt = n*nb + tsize
230 work( 1 ) = lwkopt
231 END IF
232*
233 IF( info.NE.0 ) THEN
234 CALL xerbla( 'DGEHRD', -info )
235 RETURN
236 ELSE IF( lquery ) THEN
237 RETURN
238 END IF
239*
240* Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
241*
242 DO 10 i = 1, ilo - 1
243 tau( i ) = zero
244 10 CONTINUE
245 DO 20 i = max( 1, ihi ), n - 1
246 tau( i ) = zero
247 20 CONTINUE
248*
249* Quick return if possible
250*
251 nh = ihi - ilo + 1
252 IF( nh.LE.1 ) THEN
253 work( 1 ) = 1
254 RETURN
255 END IF
256*
257* Determine the block size
258*
259 nb = min( nbmax, ilaenv( 1, 'DGEHRD', ' ', n, ilo, ihi, -1 ) )
260 nbmin = 2
261 IF( nb.GT.1 .AND. nb.LT.nh ) THEN
262*
263* Determine when to cross over from blocked to unblocked code
264* (last block is always handled by unblocked code)
265*
266 nx = max( nb, ilaenv( 3, 'DGEHRD', ' ', n, ilo, ihi, -1 ) )
267 IF( nx.LT.nh ) THEN
268*
269* Determine if workspace is large enough for blocked code
270*
271 IF( lwork.LT.n*nb+tsize ) THEN
272*
273* Not enough workspace to use optimal NB: determine the
274* minimum value of NB, and reduce NB or force use of
275* unblocked code
276*
277 nbmin = max( 2, ilaenv( 2, 'DGEHRD', ' ', n, ilo, ihi,
278 $ -1 ) )
279 IF( lwork.GE.(n*nbmin + tsize) ) THEN
280 nb = (lwork-tsize) / n
281 ELSE
282 nb = 1
283 END IF
284 END IF
285 END IF
286 END IF
287 ldwork = n
288*
289 IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
290*
291* Use unblocked code below
292*
293 i = ilo
294*
295 ELSE
296*
297* Use blocked code
298*
299 iwt = 1 + n*nb
300 DO 40 i = ilo, ihi - 1 - nx, nb
301 ib = min( nb, ihi-i )
302*
303* Reduce columns i:i+ib-1 to Hessenberg form, returning the
304* matrices V and T of the block reflector H = I - V*T*V**T
305* which performs the reduction, and also the matrix Y = A*V*T
306*
307 CALL dlahr2( ihi, i, ib, a( 1, i ), lda, tau( i ),
308 $ work( iwt ), ldt, work, ldwork )
309*
310* Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
311* right, computing A := A - Y * V**T. V(i+ib,ib-1) must be set
312* to 1
313*
314 ei = a( i+ib, i+ib-1 )
315 a( i+ib, i+ib-1 ) = one
316 CALL dgemm( 'No transpose', 'Transpose',
317 $ ihi, ihi-i-ib+1,
318 $ ib, -one, work, ldwork, a( i+ib, i ), lda, one,
319 $ a( 1, i+ib ), lda )
320 a( i+ib, i+ib-1 ) = ei
321*
322* Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
323* right
324*
325 CALL dtrmm( 'Right', 'Lower', 'Transpose',
326 $ 'Unit', i, ib-1,
327 $ one, a( i+1, i ), lda, work, ldwork )
328 DO 30 j = 0, ib-2
329 CALL daxpy( i, -one, work( ldwork*j+1 ), 1,
330 $ a( 1, i+j+1 ), 1 )
331 30 CONTINUE
332*
333* Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
334* left
335*
336 CALL dlarfb( 'Left', 'Transpose', 'Forward',
337 $ 'Columnwise',
338 $ ihi-i, n-i-ib+1, ib, a( i+1, i ), lda,
339 $ work( iwt ), ldt, a( i+1, i+ib ), lda,
340 $ work, ldwork )
341 40 CONTINUE
342 END IF
343*
344* Use unblocked code to reduce the rest of the matrix
345*
346 CALL dgehd2( n, i, ihi, a, lda, tau, work, iinfo )
347 work( 1 ) = lwkopt
348*
349 RETURN
350*
351* End of DGEHRD
352*
353 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dgehd2(n, ilo, ihi, a, lda, tau, work, info)
DGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm.
Definition dgehd2.f:149
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlahr2(n, k, nb, a, lda, tau, t, ldt, y, ldy)
DLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition dlahr2.f:181
subroutine dlarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition dlarfb.f:197
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177