LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.
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 241 of file chetrd_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 A( LDA, * ), AB( LDAB, * ),
256 $ TAU( * ), WORK( * )
257* ..
258*
259* =====================================================================
260*
261* .. Parameters ..
262 REAL RONE
263 COMPLEX ZERO, ONE, HALF
264 parameter( rone = 1.0e+0,
265 $ zero = ( 0.0e+0, 0.0e+0 ),
266 $ one = ( 1.0e+0, 0.0e+0 ),
267 $ half = ( 0.5e+0, 0.0e+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, cher2k, chemm, cgemm, ccopy,
279* ..
280* .. Intrinsic Functions ..
281 INTRINSIC min, max
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 INTEGER ILAENV2STAGE
286 REAL SROUNDUP_LWORK
288* ..
289* .. Executable Statements ..
290*
291* Determine the minimal workspace size required
292* and test the input parameters
293*
294 info = 0
295 upper = lsame( uplo, 'U' )
296 lquery = ( lwork.EQ.-1 )
297 lwmin = ilaenv2stage( 4, 'CHETRD_HE2HB', '', n, kd, -1, -1 )
298
299 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
300 info = -1
301 ELSE IF( n.LT.0 ) THEN
302 info = -2
303 ELSE IF( kd.LT.0 ) THEN
304 info = -3
305 ELSE IF( lda.LT.max( 1, n ) ) THEN
306 info = -5
307 ELSE IF( ldab.LT.max( 1, kd+1 ) ) THEN
308 info = -7
309 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
310 info = -10
311 END IF
312*
313 IF( info.NE.0 ) THEN
314 CALL xerbla( 'CHETRD_HE2HB', -info )
315 RETURN
316 ELSE IF( lquery ) THEN
317 work( 1 ) = sroundup_lwork(lwmin)
318 RETURN
319 END IF
320*
321* Quick return if possible
322* Copy the upper/lower portion of A into AB
323*
324 IF( n.LE.kd+1 ) THEN
325 IF( upper ) THEN
326 DO 100 i = 1, n
327 lk = min( kd+1, i )
328 CALL ccopy( lk, a( i-lk+1, i ), 1,
329 $ ab( kd+1-lk+1, i ), 1 )
330 100 CONTINUE
331 ELSE
332 DO 110 i = 1, n
333 lk = min( kd+1, n-i+1 )
334 CALL ccopy( lk, a( i, i ), 1, ab( 1, i ), 1 )
335 110 CONTINUE
336 ENDIF
337 work( 1 ) = 1
338 RETURN
339 END IF
340*
341* Determine the pointer position for the workspace
342*
343 ldt = kd
344 lds1 = kd
345 lt = ldt*kd
346 lw = n*kd
347 ls1 = lds1*kd
348 ls2 = lwmin - lt - lw - ls1
349* LS2 = N*MAX(KD,FACTOPTNB)
350 tpos = 1
351 wpos = tpos + lt
352 s1pos = wpos + lw
353 s2pos = s1pos + ls1
354 IF( upper ) THEN
355 ldw = kd
356 lds2 = kd
357 ELSE
358 ldw = n
359 lds2 = n
360 ENDIF
361*
362*
363* Set the workspace of the triangular matrix T to zero once such a
364* way every time T is generated the upper/lower portion will be always zero
365*
366 CALL claset( "A", ldt, kd, zero, zero, work( tpos ), ldt )
367*
368 IF( upper ) THEN
369 DO 10 i = 1, n - kd, kd
370 pn = n-i-kd+1
371 pk = min( n-i-kd+1, kd )
372*
373* Compute the LQ factorization of the current block
374*
375 CALL cgelqf( kd, pn, a( i, i+kd ), lda,
376 $ tau( i ), work( s2pos ), ls2, iinfo )
377*
378* Copy the upper portion of A into AB
379*
380 DO 20 j = i, i+pk-1
381 lk = min( kd, n-j ) + 1
382 CALL ccopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
383 20 CONTINUE
384*
385 CALL claset( 'Lower', pk, pk, zero, one,
386 $ a( i, i+kd ), lda )
387*
388* Form the matrix T
389*
390 CALL clarft( 'Forward', 'Rowwise', pn, pk,
391 $ a( i, i+kd ), lda, tau( i ),
392 $ work( tpos ), ldt )
393*
394* Compute W:
395*
396 CALL cgemm( 'Conjugate', 'No transpose', pk, pn, pk,
397 $ one, work( tpos ), ldt,
398 $ a( i, i+kd ), lda,
399 $ zero, work( s2pos ), lds2 )
400*
401 CALL chemm( 'Right', uplo, pk, pn,
402 $ one, a( i+kd, i+kd ), lda,
403 $ work( s2pos ), lds2,
404 $ zero, work( wpos ), ldw )
405*
406 CALL cgemm( 'No transpose', 'Conjugate', pk, pk, pn,
407 $ one, work( wpos ), ldw,
408 $ work( s2pos ), lds2,
409 $ zero, work( s1pos ), lds1 )
410*
411 CALL cgemm( 'No transpose', 'No transpose', pk, pn, pk,
412 $ -half, work( s1pos ), lds1,
413 $ a( i, i+kd ), lda,
414 $ one, work( wpos ), ldw )
415*
416*
417* Update the unreduced submatrix A(i+kd:n,i+kd:n), using
418* an update of the form: A := A - V'*W - W'*V
419*
420 CALL cher2k( uplo, 'Conjugate', pn, pk,
421 $ -one, a( i, i+kd ), lda,
422 $ work( wpos ), ldw,
423 $ rone, a( i+kd, i+kd ), lda )
424 10 CONTINUE
425*
426* Copy the upper band to AB which is the band storage matrix
427*
428 DO 30 j = n-kd+1, n
429 lk = min(kd, n-j) + 1
430 CALL ccopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
431 30 CONTINUE
432*
433 ELSE
434*
435* Reduce the lower triangle of A to lower band matrix
436*
437 DO 40 i = 1, n - kd, kd
438 pn = n-i-kd+1
439 pk = min( n-i-kd+1, kd )
440*
441* Compute the QR factorization of the current block
442*
443 CALL cgeqrf( pn, kd, a( i+kd, i ), lda,
444 $ tau( i ), work( s2pos ), ls2, iinfo )
445*
446* Copy the upper portion of A into AB
447*
448 DO 50 j = i, i+pk-1
449 lk = min( kd, n-j ) + 1
450 CALL ccopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
451 50 CONTINUE
452*
453 CALL claset( 'Upper', pk, pk, zero, one,
454 $ a( i+kd, i ), lda )
455*
456* Form the matrix T
457*
458 CALL clarft( 'Forward', 'Columnwise', pn, pk,
459 $ a( i+kd, i ), lda, tau( i ),
460 $ work( tpos ), ldt )
461*
462* Compute W:
463*
464 CALL cgemm( 'No transpose', 'No transpose', pn, pk, pk,
465 $ one, a( i+kd, i ), lda,
466 $ work( tpos ), ldt,
467 $ zero, work( s2pos ), lds2 )
468*
469 CALL chemm( 'Left', uplo, pn, pk,
470 $ one, a( i+kd, i+kd ), lda,
471 $ work( s2pos ), lds2,
472 $ zero, work( wpos ), ldw )
473*
474 CALL cgemm( 'Conjugate', 'No transpose', pk, pk, pn,
475 $ one, work( s2pos ), lds2,
476 $ work( wpos ), ldw,
477 $ zero, work( s1pos ), lds1 )
478*
479 CALL cgemm( 'No transpose', 'No transpose', pn, pk, pk,
480 $ -half, a( i+kd, i ), lda,
481 $ work( s1pos ), lds1,
482 $ one, work( wpos ), ldw )
483*
484*
485* Update the unreduced submatrix A(i+kd:n,i+kd:n), using
486* an update of the form: A := A - V*W' - W*V'
487*
488 CALL cher2k( uplo, 'No transpose', pn, pk,
489 $ -one, a( i+kd, i ), lda,
490 $ work( wpos ), ldw,
491 $ rone, a( i+kd, i+kd ), lda )
492* ==================================================================
493* RESTORE A FOR COMPARISON AND CHECKING TO BE REMOVED
494* DO 45 J = I, I+PK-1
495* LK = MIN( KD, N-J ) + 1
496* CALL CCOPY( LK, AB( 1, J ), 1, A( J, J ), 1 )
497* 45 CONTINUE
498* ==================================================================
499 40 CONTINUE
500*
501* Copy the lower band to AB which is the band storage matrix
502*
503 DO 60 j = n-kd+1, n
504 lk = min(kd, n-j) + 1
505 CALL ccopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
506 60 CONTINUE
507
508 END IF
509*
510 work( 1 ) = sroundup_lwork(lwmin)
511 RETURN
512*
513* End of CHETRD_HE2HB
514*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
Definition cgelqf.f:143
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
subroutine chemm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
CHEMM
Definition chemm.f:191
subroutine cher2k(uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CHER2K
Definition cher2k.f:197
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
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:163
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:106
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Here is the call graph for this function:
Here is the caller graph for this function: