LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zgehrd()

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 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 ZGEHRD
  subroutine incorporating improvements proposed by Quintana-Orti and
  Van de Geijn (2006). (See ZLAHR2.)

Definition at line 166 of file zgehrd.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*16 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*16 ZERO, ONE
186  parameter( zero = ( 0.0d+0, 0.0d+0 ),
187  $ one = ( 1.0d+0, 0.0d+0 ) )
188 * ..
189 * .. Local Scalars ..
190  LOGICAL LQUERY
191  INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
192  $ NBMIN, NH, NX
193  COMPLEX*16 EI
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL zaxpy, zgehd2, zgemm, zlahr2, zlarfb, ztrmm,
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, 'ZGEHRD', ' ', n, ilo, ihi, -1 ) )
229  lwkopt = n*nb + tsize
230  work( 1 ) = lwkopt
231  ENDIF
232 *
233  IF( info.NE.0 ) THEN
234  CALL xerbla( 'ZGEHRD', -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, 'ZGEHRD', ' ', 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, 'ZGEHRD', ' ', 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, 'ZGEHRD', ' ', 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**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 ),
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**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,
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 zgehd2( n, i, ihi, a, lda, tau, work, iinfo )
347  work( 1 ) = lwkopt
348 *
349  RETURN
350 *
351 * End of ZGEHRD
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 zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:177
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:149
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 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:181
Here is the call graph for this function:
Here is the caller graph for this function: