LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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

Download SGEHRD + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 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*
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 sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: