 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ zhetrd_he2hb()

 subroutine zhetrd_he2hb ( character UPLO, integer N, integer KD, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldab, * ) AB, integer LDAB, complex*16, dimension( * ) TAU, complex*16, dimension( * ) WORK, integer LWORK, integer INFO )

ZHETRD_HE2HB

Purpose:
``` ZHETRD_HE2HB reduces a complex Hermitian matrix A to complex Hermitian
band-diagonal form AB by a unitary similarity transformation:
Q**H * A * Q = AB.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] KD ``` KD is INTEGER The number of superdiagonals of the reduced matrix if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0. The reduced matrix is stored in the array AB.``` [in,out] A ``` A is COMPLEX*16 array, dimension (LDA,N) On entry, the Hermitian matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, if UPLO = 'U', the diagonal and first superdiagonal of A are overwritten by the corresponding elements of the tridiagonal matrix T, and the elements above the first superdiagonal, with the array TAU, represent the unitary matrix Q as a product of elementary reflectors; if UPLO = 'L', the diagonal and first subdiagonal of A are over- written by the corresponding elements of the tridiagonal matrix T, 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] AB ``` AB is COMPLEX*16 array, dimension (LDAB,N) On exit, the upper or lower triangle of the Hermitian band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).``` [in] LDAB ``` LDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.``` [out] TAU ``` TAU is COMPLEX*16 array, dimension (N-KD) The scalar factors of the elementary reflectors (see Further Details).``` [out] WORK ``` WORK is COMPLEX*16 array, dimension (LWORK) On exit, if INFO = 0, or if LWORK=-1, WORK(1) returns the size of LWORK.``` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK which should be calculated by a workspace query. LWORK = MAX(1, LWORK_QUERY) 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. LWORK_QUERY = N*KD + N*max(KD,FACTOPTNB) + 2*KD*KD where FACTOPTNB is the blocking used by the QR or LQ algorithm, usually FACTOPTNB=128 is a good choice otherwise putting LWORK=-1 will provide the size of WORK.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Further Details:
```  Implemented by Azzam Haidar.

All details are available on technical report, SC11, SC13 papers.

Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
Parallel reduction to condensed forms for symmetric eigenvalue problems
using aggregated fine-grained and memory-aware kernels. In Proceedings
of 2011 International Conference for High Performance Computing,
Networking, Storage and Analysis (SC '11), New York, NY, USA,
Article 8 , 11 pages.
http://doi.acm.org/10.1145/2063384.2063394

A. Haidar, J. Kurzak, P. Luszczek, 2013.
An improved parallel singular value algorithm and its implementation
for multicore hardware, In Proceedings of 2013 International Conference
for High Performance Computing, Networking, Storage and Analysis (SC '13).
Article 90, 12 pages.
http://doi.acm.org/10.1145/2503210.2503292

A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
A novel hybrid CPU-GPU generalized eigensolver for electronic structure
calculations based on fine-grained memory aware tasks.
International Journal of High Performance Computing Applications.
Volume 28 Issue 2, Pages 196-209, May 2014.
http://hpc.sagepub.com/content/28/2/196 ```
```  If UPLO = 'U', the matrix Q is represented as a product of elementary
reflectors

Q = H(k)**H . . . H(2)**H H(1)**H, where k = n-kd.

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+kd-1) = 0 and v(i+kd) = 1; conjg(v(i+kd+1:n)) is stored on exit in
A(i,i+kd+1:n), and tau in TAU(i).

If UPLO = 'L', the matrix Q is represented as a product of elementary
reflectors

Q = H(1) H(2) . . . H(k), where k = n-kd.

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(kd+1:i) = 0 and v(i+kd+1) = 1; v(i+kd+2:n) is stored on exit in
A(i+kd+2:n,i), and tau in TAU(i).

The contents of A on exit are illustrated by the following examples
with n = 5:

if UPLO = 'U':                       if UPLO = 'L':

(  ab  ab/v1  v1      v1     v1    )              (  ab                            )
(      ab     ab/v2   v2     v2    )              (  ab/v1  ab                     )
(             ab      ab/v3  v3    )              (  v1     ab/v2  ab              )
(                     ab     ab/v4 )              (  v1     v2     ab/v3  ab       )
(                            ab    )              (  v1     v2     v3     ab/v4 ab )

where d and e denote diagonal and off-diagonal elements of T, and vi
denotes an element of the vector defining H(i).```

Definition at line 241 of file zhetrd_he2hb.f.

243 *
244  IMPLICIT NONE
245 *
246 * -- LAPACK computational routine --
247 * -- LAPACK is a software package provided by Univ. of Tennessee, --
248 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249 *
250 * .. Scalar Arguments ..
251  CHARACTER UPLO
252  INTEGER INFO, LDA, LDAB, LWORK, N, KD
253 * ..
254 * .. Array Arguments ..
255  COMPLEX*16 A( LDA, * ), AB( LDAB, * ),
256  \$ TAU( * ), WORK( * )
257 * ..
258 *
259 * =====================================================================
260 *
261 * .. Parameters ..
262  DOUBLE PRECISION RONE
263  COMPLEX*16 ZERO, ONE, HALF
264  parameter( rone = 1.0d+0,
265  \$ zero = ( 0.0d+0, 0.0d+0 ),
266  \$ one = ( 1.0d+0, 0.0d+0 ),
267  \$ half = ( 0.5d+0, 0.0d+0 ) )
268 * ..
269 * .. Local Scalars ..
270  LOGICAL LQUERY, UPPER
271  INTEGER I, J, IINFO, LWMIN, PN, PK, LK,
272  \$ LDT, LDW, LDS2, LDS1,
273  \$ LS2, LS1, LW, LT,
274  \$ TPOS, WPOS, S2POS, S1POS
275 * ..
276 * .. External Subroutines ..
277  EXTERNAL xerbla, zher2k, zhemm, zgemm, zcopy,
279 * ..
280 * .. Intrinsic Functions ..
281  INTRINSIC min, max
282 * ..
283 * .. External Functions ..
284  LOGICAL LSAME
285  INTEGER ILAENV2STAGE
286  EXTERNAL lsame, ilaenv2stage
287 * ..
288 * .. Executable Statements ..
289 *
290 * Determine the minimal workspace size required
291 * and test the input parameters
292 *
293  info = 0
294  upper = lsame( uplo, 'U' )
295  lquery = ( lwork.EQ.-1 )
296  lwmin = ilaenv2stage( 4, 'ZHETRD_HE2HB', '', n, kd, -1, -1 )
297
298  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
299  info = -1
300  ELSE IF( n.LT.0 ) THEN
301  info = -2
302  ELSE IF( kd.LT.0 ) THEN
303  info = -3
304  ELSE IF( lda.LT.max( 1, n ) ) THEN
305  info = -5
306  ELSE IF( ldab.LT.max( 1, kd+1 ) ) THEN
307  info = -7
308  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
309  info = -10
310  END IF
311 *
312  IF( info.NE.0 ) THEN
313  CALL xerbla( 'ZHETRD_HE2HB', -info )
314  RETURN
315  ELSE IF( lquery ) THEN
316  work( 1 ) = lwmin
317  RETURN
318  END IF
319 *
320 * Quick return if possible
321 * Copy the upper/lower portion of A into AB
322 *
323  IF( n.LE.kd+1 ) THEN
324  IF( upper ) THEN
325  DO 100 i = 1, n
326  lk = min( kd+1, i )
327  CALL zcopy( lk, a( i-lk+1, i ), 1,
328  \$ ab( kd+1-lk+1, i ), 1 )
329  100 CONTINUE
330  ELSE
331  DO 110 i = 1, n
332  lk = min( kd+1, n-i+1 )
333  CALL zcopy( lk, a( i, i ), 1, ab( 1, i ), 1 )
334  110 CONTINUE
335  ENDIF
336  work( 1 ) = 1
337  RETURN
338  END IF
339 *
340 * Determine the pointer position for the workspace
341 *
342  ldt = kd
343  lds1 = kd
344  lt = ldt*kd
345  lw = n*kd
346  ls1 = lds1*kd
347  ls2 = lwmin - lt - lw - ls1
348 * LS2 = N*MAX(KD,FACTOPTNB)
349  tpos = 1
350  wpos = tpos + lt
351  s1pos = wpos + lw
352  s2pos = s1pos + ls1
353  IF( upper ) THEN
354  ldw = kd
355  lds2 = kd
356  ELSE
357  ldw = n
358  lds2 = n
359  ENDIF
360 *
361 *
362 * Set the workspace of the triangular matrix T to zero once such a
363 * way every time T is generated the upper/lower portion will be always zero
364 *
365  CALL zlaset( "A", ldt, kd, zero, zero, work( tpos ), ldt )
366 *
367  IF( upper ) THEN
368  DO 10 i = 1, n - kd, kd
369  pn = n-i-kd+1
370  pk = min( n-i-kd+1, kd )
371 *
372 * Compute the LQ factorization of the current block
373 *
374  CALL zgelqf( kd, pn, a( i, i+kd ), lda,
375  \$ tau( i ), work( s2pos ), ls2, iinfo )
376 *
377 * Copy the upper portion of A into AB
378 *
379  DO 20 j = i, i+pk-1
380  lk = min( kd, n-j ) + 1
381  CALL zcopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
382  20 CONTINUE
383 *
384  CALL zlaset( 'Lower', pk, pk, zero, one,
385  \$ a( i, i+kd ), lda )
386 *
387 * Form the matrix T
388 *
389  CALL zlarft( 'Forward', 'Rowwise', pn, pk,
390  \$ a( i, i+kd ), lda, tau( i ),
391  \$ work( tpos ), ldt )
392 *
393 * Compute W:
394 *
395  CALL zgemm( 'Conjugate', 'No transpose', pk, pn, pk,
396  \$ one, work( tpos ), ldt,
397  \$ a( i, i+kd ), lda,
398  \$ zero, work( s2pos ), lds2 )
399 *
400  CALL zhemm( 'Right', uplo, pk, pn,
401  \$ one, a( i+kd, i+kd ), lda,
402  \$ work( s2pos ), lds2,
403  \$ zero, work( wpos ), ldw )
404 *
405  CALL zgemm( 'No transpose', 'Conjugate', pk, pk, pn,
406  \$ one, work( wpos ), ldw,
407  \$ work( s2pos ), lds2,
408  \$ zero, work( s1pos ), lds1 )
409 *
410  CALL zgemm( 'No transpose', 'No transpose', pk, pn, pk,
411  \$ -half, work( s1pos ), lds1,
412  \$ a( i, i+kd ), lda,
413  \$ one, work( wpos ), ldw )
414 *
415 *
416 * Update the unreduced submatrix A(i+kd:n,i+kd:n), using
417 * an update of the form: A := A - V'*W - W'*V
418 *
419  CALL zher2k( uplo, 'Conjugate', pn, pk,
420  \$ -one, a( i, i+kd ), lda,
421  \$ work( wpos ), ldw,
422  \$ rone, a( i+kd, i+kd ), lda )
423  10 CONTINUE
424 *
425 * Copy the upper band to AB which is the band storage matrix
426 *
427  DO 30 j = n-kd+1, n
428  lk = min(kd, n-j) + 1
429  CALL zcopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
430  30 CONTINUE
431 *
432  ELSE
433 *
434 * Reduce the lower triangle of A to lower band matrix
435 *
436  DO 40 i = 1, n - kd, kd
437  pn = n-i-kd+1
438  pk = min( n-i-kd+1, kd )
439 *
440 * Compute the QR factorization of the current block
441 *
442  CALL zgeqrf( pn, kd, a( i+kd, i ), lda,
443  \$ tau( i ), work( s2pos ), ls2, iinfo )
444 *
445 * Copy the upper portion of A into AB
446 *
447  DO 50 j = i, i+pk-1
448  lk = min( kd, n-j ) + 1
449  CALL zcopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
450  50 CONTINUE
451 *
452  CALL zlaset( 'Upper', pk, pk, zero, one,
453  \$ a( i+kd, i ), lda )
454 *
455 * Form the matrix T
456 *
457  CALL zlarft( 'Forward', 'Columnwise', pn, pk,
458  \$ a( i+kd, i ), lda, tau( i ),
459  \$ work( tpos ), ldt )
460 *
461 * Compute W:
462 *
463  CALL zgemm( 'No transpose', 'No transpose', pn, pk, pk,
464  \$ one, a( i+kd, i ), lda,
465  \$ work( tpos ), ldt,
466  \$ zero, work( s2pos ), lds2 )
467 *
468  CALL zhemm( 'Left', uplo, pn, pk,
469  \$ one, a( i+kd, i+kd ), lda,
470  \$ work( s2pos ), lds2,
471  \$ zero, work( wpos ), ldw )
472 *
473  CALL zgemm( 'Conjugate', 'No transpose', pk, pk, pn,
474  \$ one, work( s2pos ), lds2,
475  \$ work( wpos ), ldw,
476  \$ zero, work( s1pos ), lds1 )
477 *
478  CALL zgemm( 'No transpose', 'No transpose', pn, pk, pk,
479  \$ -half, a( i+kd, i ), lda,
480  \$ work( s1pos ), lds1,
481  \$ one, work( wpos ), ldw )
482 *
483 *
484 * Update the unreduced submatrix A(i+kd:n,i+kd:n), using
485 * an update of the form: A := A - V*W' - W*V'
486 *
487  CALL zher2k( uplo, 'No transpose', pn, pk,
488  \$ -one, a( i+kd, i ), lda,
489  \$ work( wpos ), ldw,
490  \$ rone, a( i+kd, i+kd ), lda )
491 * ==================================================================
492 * RESTORE A FOR COMPARISON AND CHECKING TO BE REMOVED
493 * DO 45 J = I, I+PK-1
494 * LK = MIN( KD, N-J ) + 1
495 * CALL ZCOPY( LK, AB( 1, J ), 1, A( J, J ), 1 )
496 * 45 CONTINUE
497 * ==================================================================
498  40 CONTINUE
499 *
500 * Copy the lower band to AB which is the band storage matrix
501 *
502  DO 60 j = n-kd+1, n
503  lk = min(kd, n-j) + 1
504  CALL zcopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
505  60 CONTINUE
506
507  END IF
508 *
509  work( 1 ) = lwmin
510  RETURN
511 *
512 * End of ZHETRD_HE2HB
513 *
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:149
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zhemm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZHEMM
Definition: zhemm.f:191
subroutine zher2k(UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZHER2K
Definition: zher2k.f:198
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:143
subroutine zlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition: zlarft.f:163
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
Here is the call graph for this function:
Here is the caller graph for this function: