LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ chetrd_he2hb()

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

CHETRD_HE2HB

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

Purpose:
 CHETRD_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 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 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 array, dimension (N-KD)
          The scalar factors of the elementary reflectors (see Further
          Details).
[out]WORK
          WORK is COMPLEX 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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2017
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).
  Denver, Colorado, USA, 2013.
  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 245 of file chetrd_he2hb.f.

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