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

◆ cgehrd()

subroutine cgehrd ( integer  n,
integer  ilo,
integer  ihi,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( * )  tau,
complex, dimension( * )  work,
integer  lwork,
integer  info 
)

CGEHRD

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

Purpose:
 CGEHRD reduces a complex general matrix A to upper Hessenberg form H by
 an unitary similarity transformation:  Q**H * 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 CGEBAL; 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 COMPLEX 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 unitary 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 COMPLEX 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 COMPLEX 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**H

  where tau is a complex scalar, and v is a complex 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 CGEHRD
  subroutine incorporating improvements proposed by Quintana-Orti and
  Van de Geijn (2006). (See CLAHR2.)

Definition at line 166 of file cgehrd.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 COMPLEX 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 COMPLEX ZERO, ONE
186 parameter( zero = ( 0.0e+0, 0.0e+0 ),
187 $ one = ( 1.0e+0, 0.0e+0 ) )
188* ..
189* .. Local Scalars ..
190 LOGICAL LQUERY
191 INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
192 $ NBMIN, NH, NX
193 COMPLEX EI
194* ..
195* .. External Subroutines ..
196 EXTERNAL caxpy, cgehd2, cgemm, clahr2, clarfb, ctrmm,
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, 'CGEHRD', ' ', 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( 'CGEHRD', -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, 'CGEHRD', ' ', 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, 'CGEHRD', ' ', 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, 'CGEHRD', ' ', 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**H
306* which performs the reduction, and also the matrix Y = A*V*T
307*
308 CALL clahr2( 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**H. 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 cgemm( 'No transpose', 'Conjugate 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 ctrmm( 'Right', 'Lower', 'Conjugate transpose',
327 $ 'Unit', i, ib-1,
328 $ one, a( i+1, i ), lda, work, ldwork )
329 DO 30 j = 0, ib-2
330 CALL caxpy( 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 clarfb( 'Left', 'Conjugate 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 cgehd2( n, i, ihi, a, lda, tau, work, iinfo )
348 work( 1 ) = sroundup_lwork(lwkopt)
349*
350 RETURN
351*
352* End of CGEHRD
353*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine cgehd2(n, ilo, ihi, a, lda, tau, work, info)
CGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm.
Definition cgehd2.f:149
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine clahr2(n, k, nb, a, lda, tau, t, ldt, y, ldy)
CLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition clahr2.f:181
subroutine clarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
Definition clarfb.f:197
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: