LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
zgehrd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine zgehrd (N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
 ZGEHRD More...
 

Function/Subroutine Documentation

subroutine zgehrd ( integer  N,
integer  ILO,
integer  IHI,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZGEHRD

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

Purpose:
 ZGEHRD 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 ZGEBAL; 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*16 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*16 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*16 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 optimum performance LWORK >= N*NB, where NB is the
          optimal blocksize.

          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 2011
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 170 of file zgehrd.f.

170 *
171 * -- LAPACK computational routine (version 3.4.0) --
172 * -- LAPACK is a software package provided by Univ. of Tennessee, --
173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174 * November 2011
175 *
176 * .. Scalar Arguments ..
177  INTEGER ihi, ilo, info, lda, lwork, n
178 * ..
179 * .. Array Arguments ..
180  COMPLEX*16 a( lda, * ), tau( * ), work( * )
181 * ..
182 *
183 * =====================================================================
184 *
185 * .. Parameters ..
186  INTEGER nbmax, ldt
187  parameter( nbmax = 64, ldt = nbmax+1 )
188  COMPLEX*16 zero, one
189  parameter( zero = ( 0.0d+0, 0.0d+0 ),
190  $ one = ( 1.0d+0, 0.0d+0 ) )
191 * ..
192 * .. Local Scalars ..
193  LOGICAL lquery
194  INTEGER i, ib, iinfo, iws, j, ldwork, lwkopt, nb,
195  $ nbmin, nh, nx
196  COMPLEX*16 ei
197 * ..
198 * .. Local Arrays ..
199  COMPLEX*16 t( ldt, nbmax )
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL zaxpy, zgehd2, zgemm, zlahr2, zlarfb, ztrmm,
203  $ xerbla
204 * ..
205 * .. Intrinsic Functions ..
206  INTRINSIC max, min
207 * ..
208 * .. External Functions ..
209  INTEGER ilaenv
210  EXTERNAL ilaenv
211 * ..
212 * .. Executable Statements ..
213 *
214 * Test the input parameters
215 *
216  info = 0
217  nb = min( nbmax, ilaenv( 1, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
218  lwkopt = n*nb
219  work( 1 ) = lwkopt
220  lquery = ( lwork.EQ.-1 )
221  IF( n.LT.0 ) THEN
222  info = -1
223  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
224  info = -2
225  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
226  info = -3
227  ELSE IF( lda.LT.max( 1, n ) ) THEN
228  info = -5
229  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
230  info = -8
231  END IF
232  IF( info.NE.0 ) THEN
233  CALL xerbla( 'ZGEHRD', -info )
234  RETURN
235  ELSE IF( lquery ) THEN
236  RETURN
237  END IF
238 *
239 * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
240 *
241  DO 10 i = 1, ilo - 1
242  tau( i ) = zero
243  10 CONTINUE
244  DO 20 i = max( 1, ihi ), n - 1
245  tau( i ) = zero
246  20 CONTINUE
247 *
248 * Quick return if possible
249 *
250  nh = ihi - ilo + 1
251  IF( nh.LE.1 ) THEN
252  work( 1 ) = 1
253  RETURN
254  END IF
255 *
256 * Determine the block size
257 *
258  nb = min( nbmax, ilaenv( 1, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
259  nbmin = 2
260  iws = 1
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, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
267  IF( nx.LT.nh ) THEN
268 *
269 * Determine if workspace is large enough for blocked code
270 *
271  iws = n*nb
272  IF( lwork.LT.iws ) 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, 'ZGEHRD', ' ', n, ilo, ihi,
279  $ -1 ) )
280  IF( lwork.GE.n*nbmin ) THEN
281  nb = lwork / 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  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**H
305 * which performs the reduction, and also the matrix Y = A*V*T
306 *
307  CALL zlahr2( ihi, i, ib, a( 1, i ), lda, tau( i ), t, ldt,
308  $ 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**H. 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 zgemm( 'No transpose', 'Conjugate 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 ztrmm( 'Right', 'Lower', 'Conjugate transpose',
326  $ 'Unit', i, ib-1,
327  $ one, a( i+1, i ), lda, work, ldwork )
328  DO 30 j = 0, ib-2
329  CALL zaxpy( 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 zlarfb( 'Left', 'Conjugate transpose', 'Forward',
337  $ 'Columnwise',
338  $ ihi-i, n-i-ib+1, ib, a( i+1, i ), lda, t, ldt,
339  $ a( i+1, i+ib ), lda, work, ldwork )
340  40 CONTINUE
341  END IF
342 *
343 * Use unblocked code to reduce the rest of the matrix
344 *
345  CALL zgehd2( n, i, ihi, a, lda, tau, work, iinfo )
346  work( 1 ) = iws
347 *
348  RETURN
349 *
350 * End of ZGEHRD
351 *
subroutine zlahr2(N, K, NB, A, LDA, TAU, T, LDT, Y, LDY)
ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elemen...
Definition: zlahr2.f:183
subroutine zlarfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix...
Definition: zlarfb.f:197
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine zgehd2(N, ILO, IHI, A, LDA, TAU, WORK, INFO)
ZGEHD2 reduces a general square matrix to upper Hessenberg form using an unblocked algorithm...
Definition: zgehd2.f:151
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:179

Here is the call graph for this function:

Here is the caller graph for this function: