LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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```
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 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 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, 'CHETRD_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( 'CHETRD_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 ccopy( 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 ccopy( 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 claset( "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 cgelqf( 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 ccopy( lk, a( j, j ), lda, ab( kd+1, j ), ldab-1 )
382 20 CONTINUE
383*
384 CALL claset( 'Lower', pk, pk, zero, one,
385 \$ a( i, i+kd ), lda )
386*
387* Form the matrix T
388*
389 CALL clarft( 'Forward', 'Rowwise', pn, pk,
390 \$ a( i, i+kd ), lda, tau( i ),
391 \$ work( tpos ), ldt )
392*
393* Compute W:
394*
395 CALL cgemm( '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 chemm( 'Right', uplo, pk, pn,
401 \$ one, a( i+kd, i+kd ), lda,
402 \$ work( s2pos ), lds2,
403 \$ zero, work( wpos ), ldw )
404*
405 CALL cgemm( 'No transpose', 'Conjugate', pk, pk, pn,
406 \$ one, work( wpos ), ldw,
407 \$ work( s2pos ), lds2,
408 \$ zero, work( s1pos ), lds1 )
409*
410 CALL cgemm( '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 cher2k( 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 ccopy( 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 cgeqrf( 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 ccopy( lk, a( j, j ), 1, ab( 1, j ), 1 )
450 50 CONTINUE
451*
452 CALL claset( 'Upper', pk, pk, zero, one,
453 \$ a( i+kd, i ), lda )
454*
455* Form the matrix T
456*
457 CALL clarft( 'Forward', 'Columnwise', pn, pk,
458 \$ a( i+kd, i ), lda, tau( i ),
459 \$ work( tpos ), ldt )
460*
461* Compute W:
462*
463 CALL cgemm( '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 chemm( 'Left', uplo, pn, pk,
469 \$ one, a( i+kd, i+kd ), lda,
470 \$ work( s2pos ), lds2,
471 \$ zero, work( wpos ), ldw )
472*
473 CALL cgemm( 'Conjugate', 'No transpose', pk, pk, pn,
474 \$ one, work( s2pos ), lds2,
475 \$ work( wpos ), ldw,
476 \$ zero, work( s1pos ), lds1 )
477*
478 CALL cgemm( '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 cher2k( 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 CCOPY( 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 ccopy( 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 CHETRD_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 ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
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
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:146
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:143
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
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
Here is the call graph for this function:
Here is the caller graph for this function: