LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ sgehrd()

 subroutine sgehrd ( integer N, integer ILO, integer IHI, real, dimension( lda, * ) A, integer LDA, real, dimension( * ) TAU, real, dimension( * ) WORK, integer LWORK, integer INFO )

SGEHRD

Purpose:
SGEHRD reduces a real general matrix A to upper Hessenberg form H by
an orthogonal similarity transformation:  Q**T * A * Q = H .
Parameters
 [in] N N is INTEGER The order of the matrix A. N >= 0. [in] ILO ILO is INTEGER [in] IHI IHI is INTEGER It is assumed that A is already upper triangular in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally set by a previous call to SGEBAL; otherwise they should be set to 1 and N respectively. See Further Details. 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. [in,out] A A is REAL array, dimension (LDA,N) On entry, the N-by-N general matrix to be reduced. On exit, the upper triangle and the first subdiagonal of A are overwritten with the upper Hessenberg matrix H, and the elements below the first subdiagonal, with the array TAU, represent the orthogonal matrix Q as a product of elementary reflectors. See Further Details. [in] LDA LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N). [out] TAU TAU is REAL array, dimension (N-1) The scalar factors of the elementary reflectors (see Further Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to zero. [out] WORK WORK is REAL array, dimension (LWORK) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. [in] LWORK LWORK is INTEGER The length of the array WORK. LWORK >= max(1,N). For good performance, LWORK should generally be larger. If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA. [out] INFO INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value.
Further Details:
The matrix Q is represented as a product of (ihi-ilo) elementary
reflectors

Q = H(ilo) H(ilo+1) . . . H(ihi-1).

Each H(i) has the form

H(i) = I - tau * v * v**T

where tau is a real scalar, and v is a real vector with
v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
exit in A(i+2:ihi,i), and tau in TAU(i).

The contents of A are illustrated by the following example, with
n = 7, ilo = 2 and ihi = 6:

on entry,                        on exit,

( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
(     a   a   a   a   a   a )    (      a   h   h   h   h   a )
(     a   a   a   a   a   a )    (      h   h   h   h   h   h )
(     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
(     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
(     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
(                         a )    (                          a )

where a denotes an element of the original matrix A, h denotes a
modified element of the upper Hessenberg matrix H, and vi denotes an
element of the vector defining H(i).

This file is a slight modification of LAPACK-3.0's SGEHRD
subroutine incorporating improvements proposed by Quintana-Orti and
Van de Geijn (2006). (See SLAHR2.)

Definition at line 166 of file sgehrd.f.

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  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, 'SGEHRD', ' ', 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( 'SGEHRD', -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, 'SGEHRD', ' ', 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, 'SGEHRD', ' ', 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, 'SGEHRD', ' ', 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 slahr2( 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 sgemm( '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 strmm( '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 saxpy( 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 slarfb( '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 sgehd2( n, i, ihi, a, lda, tau, work, iinfo )
347  work( 1 ) = lwkopt
348 *
349  RETURN
350 *
351 * End of SGEHRD
352 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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 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 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 saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:177
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: