LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgehrd.f
Go to the documentation of this file.
1*> \brief \b SGEHRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGEHRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgehrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgehrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgehrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGEHRD( 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* REAL A( LDA, * ), TAU( * ), WORK( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> SGEHRD 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 SGEBAL; 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 REAL 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 REAL 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 REAL 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 SGEHRD
161*> subroutine incorporating improvements proposed by Quintana-Orti and
162*> Van de Geijn (2006). (See SLAHR2.)
163*> \endverbatim
164*>
165* =====================================================================
166 SUBROUTINE sgehrd( 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 REAL 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 REAL ZERO, ONE
186 parameter( zero = 0.0e+0,
187 $ one = 1.0e+0 )
188* ..
189* .. Local Scalars ..
190 LOGICAL LQUERY
191 INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
192 $ NBMIN, NH, NX
193 REAL EI
194* ..
195* .. External Subroutines ..
196 EXTERNAL saxpy, sgehd2, sgemm, slahr2, slarfb, strmm,
197 $ xerbla
198* ..
199* .. Intrinsic Functions ..
200 INTRINSIC max, min
201* ..
202* .. External Functions ..
203 INTEGER ILAENV
204 REAL SROUNDUP_LWORK
205 EXTERNAL ilaenv, sroundup_lwork
206* ..
207* .. Executable Statements ..
208*
209* Test the input parameters
210*
211 info = 0
212 lquery = ( lwork.EQ.-1 )
213 IF( n.LT.0 ) THEN
214 info = -1
215 ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
216 info = -2
217 ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
218 info = -3
219 ELSE IF( lda.LT.max( 1, n ) ) THEN
220 info = -5
221 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
222 info = -8
223 END IF
224*
225 IF( info.EQ.0 ) THEN
226*
227* Compute the workspace requirements
228*
229 nb = min( nbmax, ilaenv( 1, 'SGEHRD', ' ', n, ilo, ihi, -1 ) )
230 lwkopt = n*nb + tsize
231 work( 1 ) = sroundup_lwork(lwkopt)
232 END IF
233*
234 IF( info.NE.0 ) THEN
235 CALL xerbla( 'SGEHRD', -info )
236 RETURN
237 ELSE IF( lquery ) THEN
238 RETURN
239 END IF
240*
241* Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
242*
243 DO 10 i = 1, ilo - 1
244 tau( i ) = zero
245 10 CONTINUE
246 DO 20 i = max( 1, ihi ), n - 1
247 tau( i ) = zero
248 20 CONTINUE
249*
250* Quick return if possible
251*
252 nh = ihi - ilo + 1
253 IF( nh.LE.1 ) THEN
254 work( 1 ) = 1
255 RETURN
256 END IF
257*
258* Determine the block size
259*
260 nb = min( nbmax, ilaenv( 1, 'SGEHRD', ' ', n, ilo, ihi, -1 ) )
261 nbmin = 2
262 IF( nb.GT.1 .AND. nb.LT.nh ) THEN
263*
264* Determine when to cross over from blocked to unblocked code
265* (last block is always handled by unblocked code)
266*
267 nx = max( nb, ilaenv( 3, 'SGEHRD', ' ', n, ilo, ihi, -1 ) )
268 IF( nx.LT.nh ) THEN
269*
270* Determine if workspace is large enough for blocked code
271*
272 IF( lwork.LT.n*nb+tsize ) THEN
273*
274* Not enough workspace to use optimal NB: determine the
275* minimum value of NB, and reduce NB or force use of
276* unblocked code
277*
278 nbmin = max( 2, ilaenv( 2, 'SGEHRD', ' ', n, ilo, ihi,
279 $ -1 ) )
280 IF( lwork.GE.(n*nbmin + tsize) ) THEN
281 nb = (lwork-tsize) / n
282 ELSE
283 nb = 1
284 END IF
285 END IF
286 END IF
287 END IF
288 ldwork = n
289*
290 IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
291*
292* Use unblocked code below
293*
294 i = ilo
295*
296 ELSE
297*
298* Use blocked code
299*
300 iwt = 1 + n*nb
301 DO 40 i = ilo, ihi - 1 - nx, nb
302 ib = min( nb, ihi-i )
303*
304* Reduce columns i:i+ib-1 to Hessenberg form, returning the
305* matrices V and T of the block reflector H = I - V*T*V**T
306* which performs the reduction, and also the matrix Y = A*V*T
307*
308 CALL slahr2( ihi, i, ib, a( 1, i ), lda, tau( i ),
309 $ work( iwt ), ldt, work, ldwork )
310*
311* Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
312* right, computing A := A - Y * V**T. V(i+ib,ib-1) must be set
313* to 1
314*
315 ei = a( i+ib, i+ib-1 )
316 a( i+ib, i+ib-1 ) = one
317 CALL sgemm( 'No transpose', 'Transpose',
318 $ ihi, ihi-i-ib+1,
319 $ ib, -one, work, ldwork, a( i+ib, i ), lda, one,
320 $ a( 1, i+ib ), lda )
321 a( i+ib, i+ib-1 ) = ei
322*
323* Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
324* right
325*
326 CALL strmm( 'Right', 'Lower', 'Transpose',
327 $ 'Unit', i, ib-1,
328 $ one, a( i+1, i ), lda, work, ldwork )
329 DO 30 j = 0, ib-2
330 CALL saxpy( i, -one, work( ldwork*j+1 ), 1,
331 $ a( 1, i+j+1 ), 1 )
332 30 CONTINUE
333*
334* Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
335* left
336*
337 CALL slarfb( 'Left', 'Transpose', 'Forward',
338 $ 'Columnwise',
339 $ ihi-i, n-i-ib+1, ib, a( i+1, i ), lda,
340 $ work( iwt ), ldt, a( i+1, i+ib ), lda,
341 $ work, ldwork )
342 40 CONTINUE
343 END IF
344*
345* Use unblocked code to reduce the rest of the matrix
346*
347 CALL sgehd2( n, i, ihi, a, lda, tau, work, iinfo )
348 work( 1 ) = sroundup_lwork(lwkopt)
349*
350 RETURN
351*
352* End of SGEHRD
353*
354 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
subroutine sgehd2(n, ilo, ihi, a, lda, tau, work, info)
SGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm.
Definition sgehd2.f:149
subroutine sgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
SGEHRD
Definition sgehrd.f:167
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slahr2(n, k, nb, a, lda, tau, t, ldt, y, ldy)
SLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition slahr2.f:181
subroutine slarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition slarfb.f:197
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177