LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2015
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 DGEHRD
  subroutine incorporating improvements proposed by Quintana-Orti and
  Van de Geijn (2006). (See DLAHR2.)

Definition at line 169 of file cgehrd.f.

169 *
170 * -- LAPACK computational routine (version 3.6.0) --
171 * -- LAPACK is a software package provided by Univ. of Tennessee, --
172 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173 * November 2015
174 *
175 * .. Scalar Arguments ..
176  INTEGER ihi, ilo, info, lda, lwork, n
177 * ..
178 * .. Array Arguments ..
179  COMPLEX a( lda, * ), tau( * ), work( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  INTEGER nbmax, ldt, tsize
186  parameter ( nbmax = 64, ldt = nbmax+1,
187  $ tsize = ldt*nbmax )
188  COMPLEX zero, one
189  parameter ( zero = ( 0.0e+0, 0.0e+0 ),
190  $ one = ( 1.0e+0, 0.0e+0 ) )
191 * ..
192 * .. Local Scalars ..
193  LOGICAL lquery
194  INTEGER i, ib, iinfo, iwt, j, ldwork, lwkopt, nb,
195  $ nbmin, nh, nx
196  COMPLEX ei
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL caxpy, cgehd2, cgemm, clahr2, clarfb, ctrmm,
200  $ xerbla
201 * ..
202 * .. Intrinsic Functions ..
203  INTRINSIC max, min
204 * ..
205 * .. External Functions ..
206  INTEGER ilaenv
207  EXTERNAL ilaenv
208 * ..
209 * .. Executable Statements ..
210 *
211 * Test the input parameters
212 *
213  info = 0
214  lquery = ( lwork.EQ.-1 )
215  IF( n.LT.0 ) THEN
216  info = -1
217  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
218  info = -2
219  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
220  info = -3
221  ELSE IF( lda.LT.max( 1, n ) ) THEN
222  info = -5
223  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
224  info = -8
225  END IF
226 *
227  IF( info.EQ.0 ) THEN
228 *
229 * Compute the workspace requirements
230 *
231  nb = min( nbmax, ilaenv( 1, 'CGEHRD', ' ', n, ilo, ihi, -1 ) )
232  lwkopt = n*nb + tsize
233  work( 1 ) = lwkopt
234  END IF
235 *
236  IF( info.NE.0 ) THEN
237  CALL xerbla( 'CGEHRD', -info )
238  RETURN
239  ELSE IF( lquery ) THEN
240  RETURN
241  END IF
242 *
243 * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
244 *
245  DO 10 i = 1, ilo - 1
246  tau( i ) = zero
247  10 CONTINUE
248  DO 20 i = max( 1, ihi ), n - 1
249  tau( i ) = zero
250  20 CONTINUE
251 *
252 * Quick return if possible
253 *
254  nh = ihi - ilo + 1
255  IF( nh.LE.1 ) THEN
256  work( 1 ) = 1
257  RETURN
258  END IF
259 *
260 * Determine the block size
261 *
262  nb = min( nbmax, ilaenv( 1, 'CGEHRD', ' ', n, ilo, ihi, -1 ) )
263  nbmin = 2
264  IF( nb.GT.1 .AND. nb.LT.nh ) THEN
265 *
266 * Determine when to cross over from blocked to unblocked code
267 * (last block is always handled by unblocked code)
268 *
269  nx = max( nb, ilaenv( 3, 'CGEHRD', ' ', n, ilo, ihi, -1 ) )
270  IF( nx.LT.nh ) THEN
271 *
272 * Determine if workspace is large enough for blocked code
273 *
274  IF( lwork.LT.n*nb+tsize ) THEN
275 *
276 * Not enough workspace to use optimal NB: determine the
277 * minimum value of NB, and reduce NB or force use of
278 * unblocked code
279 *
280  nbmin = max( 2, ilaenv( 2, 'CGEHRD', ' ', n, ilo, ihi,
281  $ -1 ) )
282  IF( lwork.GE.(n*nbmin+tsize) ) THEN
283  nb = (lwork-tsize) / n
284  ELSE
285  nb = 1
286  END IF
287  END IF
288  END IF
289  END IF
290  ldwork = n
291 *
292  IF( nb.LT.nbmin .OR. nb.GE.nh ) THEN
293 *
294 * Use unblocked code below
295 *
296  i = ilo
297 *
298  ELSE
299 *
300 * Use blocked code
301 *
302  iwt = 1 + n*nb
303  DO 40 i = ilo, ihi - 1 - nx, nb
304  ib = min( nb, ihi-i )
305 *
306 * Reduce columns i:i+ib-1 to Hessenberg form, returning the
307 * matrices V and T of the block reflector H = I - V*T*V**H
308 * which performs the reduction, and also the matrix Y = A*V*T
309 *
310  CALL clahr2( ihi, i, ib, a( 1, i ), lda, tau( i ),
311  $ work( iwt ), ldt, work, ldwork )
312 *
313 * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
314 * right, computing A := A - Y * V**H. V(i+ib,ib-1) must be set
315 * to 1
316 *
317  ei = a( i+ib, i+ib-1 )
318  a( i+ib, i+ib-1 ) = one
319  CALL cgemm( 'No transpose', 'Conjugate transpose',
320  $ ihi, ihi-i-ib+1,
321  $ ib, -one, work, ldwork, a( i+ib, i ), lda, one,
322  $ a( 1, i+ib ), lda )
323  a( i+ib, i+ib-1 ) = ei
324 *
325 * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
326 * right
327 *
328  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
329  $ 'Unit', i, ib-1,
330  $ one, a( i+1, i ), lda, work, ldwork )
331  DO 30 j = 0, ib-2
332  CALL caxpy( i, -one, work( ldwork*j+1 ), 1,
333  $ a( 1, i+j+1 ), 1 )
334  30 CONTINUE
335 *
336 * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
337 * left
338 *
339  CALL clarfb( 'Left', 'Conjugate transpose', 'Forward',
340  $ 'Columnwise',
341  $ ihi-i, n-i-ib+1, ib, a( i+1, i ), lda,
342  $ work( iwt ), ldt, a( i+1, i+ib ), lda,
343  $ work, ldwork )
344  40 CONTINUE
345  END IF
346 *
347 * Use unblocked code to reduce the rest of the matrix
348 *
349  CALL cgehd2( n, i, ihi, a, lda, tau, work, iinfo )
350  work( 1 ) = lwkopt
351 *
352  RETURN
353 *
354 * End of CGEHRD
355 *
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
clahr2
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:183
cgehd2
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:151
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
ctrmm
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
clarfb
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
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189

Here is the call graph for this function:

Here is the caller graph for this function: