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

◆ chetrd_hb2st()

subroutine chetrd_hb2st ( character  stage1,
character  vect,
character  uplo,
integer  n,
integer  kd,
complex, dimension( ldab, * )  ab,
integer  ldab,
real, dimension( * )  d,
real, dimension( * )  e,
complex, dimension( * )  hous,
integer  lhous,
complex, dimension( * )  work,
integer  lwork,
integer  info 
)

CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T

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

Purpose:
 CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric
 tridiagonal form T by a unitary similarity transformation:
 Q**H * A * Q = T.
Parameters
[in]STAGE1
          STAGE1 is CHARACTER*1
          = 'N':  "No": to mention that the stage 1 of the reduction
                  from dense to band using the chetrd_he2hb routine
                  was not called before this routine to reproduce AB.
                  In other term this routine is called as standalone.
          = 'Y':  "Yes": to mention that the stage 1 of the
                  reduction from dense to band using the chetrd_he2hb
                  routine has been called to produce AB (e.g., AB is
                  the output of chetrd_he2hb.
[in]VECT
          VECT is CHARACTER*1
          = 'N':  No need for the Housholder representation,
                  and thus LHOUS is of size max(1, 4*N);
          = 'V':  the Householder representation is needed to
                  either generate or to apply Q later on,
                  then LHOUS is to be queried and computed.
                  (NOT AVAILABLE IN THIS RELEASE).
[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 matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
[in,out]AB
          AB is COMPLEX array, dimension (LDAB,N)
          On entry, 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).
          On exit, the diagonal elements of AB are overwritten by the
          diagonal elements of the tridiagonal matrix T; if KD > 0, the
          elements on the first superdiagonal (if UPLO = 'U') or the
          first subdiagonal (if UPLO = 'L') are overwritten by the
          off-diagonal elements of T; the rest of AB is overwritten by
          values generated during the reduction.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD+1.
[out]D
          D is REAL array, dimension (N)
          The diagonal elements of the tridiagonal matrix T.
[out]E
          E is REAL array, dimension (N-1)
          The off-diagonal elements of the tridiagonal matrix T:
          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
[out]HOUS
          HOUS is COMPLEX array, dimension LHOUS, that
          store the Householder representation.
[in]LHOUS
          LHOUS is INTEGER
          The dimension of the array HOUS. LHOUS = MAX(1, dimension)
          If LWORK = -1, or LHOUS=-1,
          then a query is assumed; the routine
          only calculates the optimal size of the HOUS array, returns
          this value as the first entry of the HOUS array, and no error
          message related to LHOUS is issued by XERBLA.
          LHOUS = MAX(1, dimension) where
          dimension = 4*N if VECT='N'
          not available now if VECT='H'
[out]WORK
          WORK is COMPLEX array, dimension LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK = MAX(1, dimension)
          If LWORK = -1, or LHOUS=-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 = MAX(1, dimension) where
          dimension   = (2KD+1)*N + KD*NTHREADS
          where KD is the blocking size of the reduction,
          FACTOPTNB is the blocking used by the QR or LQ
          algorithm, usually FACTOPTNB=128 is a good choice
          NTHREADS is the number of threads used when
          openMP compilation is enabled, otherwise =1.
[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

Definition at line 228 of file chetrd_hb2st.F.

230*
231*
232#if defined(_OPENMP)
233 use omp_lib
234#endif
235*
236 IMPLICIT NONE
237*
238* -- LAPACK computational routine --
239* -- LAPACK is a software package provided by Univ. of Tennessee, --
240* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
241*
242* .. Scalar Arguments ..
243 CHARACTER STAGE1, UPLO, VECT
244 INTEGER N, KD, LDAB, LHOUS, LWORK, INFO
245* ..
246* .. Array Arguments ..
247 REAL D( * ), E( * )
248 COMPLEX AB( LDAB, * ), HOUS( * ), WORK( * )
249* ..
250*
251* =====================================================================
252*
253* .. Parameters ..
254 REAL RZERO
255 COMPLEX ZERO, ONE
256 parameter( rzero = 0.0e+0,
257 $ zero = ( 0.0e+0, 0.0e+0 ),
258 $ one = ( 1.0e+0, 0.0e+0 ) )
259* ..
260* .. Local Scalars ..
261 LOGICAL LQUERY, WANTQ, UPPER, AFTERS1
262 INTEGER I, M, K, IB, SWEEPID, MYID, SHIFT, STT, ST,
263 $ ED, STIND, EDIND, BLKLASTIND, COLPT, THED,
264 $ STEPERCOL, GRSIZ, THGRSIZ, THGRNB, THGRID,
265 $ NBTILES, TTYPE, TID, NTHREADS, DEBUG,
266 $ ABDPOS, ABOFDPOS, DPOS, OFDPOS, AWPOS,
267 $ INDA, INDW, APOS, SIZEA, LDA, INDV, INDTAU,
268 $ SICEV, SIZETAU, LDV, LHMIN, LWMIN
269 REAL ABSTMP
270 COMPLEX TMP
271* ..
272* .. External Subroutines ..
274* ..
275* .. Intrinsic Functions ..
276 INTRINSIC min, max, ceiling, real
277* ..
278* .. External Functions ..
279 LOGICAL LSAME
280 INTEGER ILAENV2STAGE
281 REAL SROUNDUP_LWORK
283* ..
284* .. Executable Statements ..
285*
286* Determine the minimal workspace size required.
287* Test the input parameters
288*
289 debug = 0
290 info = 0
291 afters1 = lsame( stage1, 'Y' )
292 wantq = lsame( vect, 'V' )
293 upper = lsame( uplo, 'U' )
294 lquery = ( lwork.EQ.-1 ) .OR. ( lhous.EQ.-1 )
295*
296* Determine the block size, the workspace size and the hous size.
297*
298 ib = ilaenv2stage( 2, 'CHETRD_HB2ST', vect, n, kd, -1, -1 )
299 lhmin = ilaenv2stage( 3, 'CHETRD_HB2ST', vect, n, kd, ib, -1 )
300 lwmin = ilaenv2stage( 4, 'CHETRD_HB2ST', vect, n, kd, ib, -1 )
301*
302 IF( .NOT.afters1 .AND. .NOT.lsame( stage1, 'N' ) ) THEN
303 info = -1
304 ELSE IF( .NOT.lsame( vect, 'N' ) ) THEN
305 info = -2
306 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
307 info = -3
308 ELSE IF( n.LT.0 ) THEN
309 info = -4
310 ELSE IF( kd.LT.0 ) THEN
311 info = -5
312 ELSE IF( ldab.LT.(kd+1) ) THEN
313 info = -7
314 ELSE IF( lhous.LT.lhmin .AND. .NOT.lquery ) THEN
315 info = -11
316 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
317 info = -13
318 END IF
319*
320 IF( info.EQ.0 ) THEN
321 hous( 1 ) = lhmin
322 work( 1 ) = sroundup_lwork(lwmin)
323 END IF
324*
325 IF( info.NE.0 ) THEN
326 CALL xerbla( 'CHETRD_HB2ST', -info )
327 RETURN
328 ELSE IF( lquery ) THEN
329 RETURN
330 END IF
331*
332* Quick return if possible
333*
334 IF( n.EQ.0 ) THEN
335 hous( 1 ) = 1
336 work( 1 ) = 1
337 RETURN
338 END IF
339*
340* Determine pointer position
341*
342 ldv = kd + ib
343 sizetau = 2 * n
344 sicev = 2 * n
345 indtau = 1
346 indv = indtau + sizetau
347 lda = 2 * kd + 1
348 sizea = lda * n
349 inda = 1
350 indw = inda + sizea
351 nthreads = 1
352 tid = 0
353*
354 IF( upper ) THEN
355 apos = inda + kd
356 awpos = inda
357 dpos = apos + kd
358 ofdpos = dpos - 1
359 abdpos = kd + 1
360 abofdpos = kd
361 ELSE
362 apos = inda
363 awpos = inda + kd + 1
364 dpos = apos
365 ofdpos = dpos + 1
366 abdpos = 1
367 abofdpos = 2
368
369 ENDIF
370*
371* Case KD=0:
372* The matrix is diagonal. We just copy it (convert to "real" for
373* complex because D is double and the imaginary part should be 0)
374* and store it in D. A sequential code here is better or
375* in a parallel environment it might need two cores for D and E
376*
377 IF( kd.EQ.0 ) THEN
378 DO 30 i = 1, n
379 d( i ) = real( ab( abdpos, i ) )
380 30 CONTINUE
381 DO 40 i = 1, n-1
382 e( i ) = rzero
383 40 CONTINUE
384*
385 hous( 1 ) = 1
386 work( 1 ) = 1
387 RETURN
388 END IF
389*
390* Case KD=1:
391* The matrix is already Tridiagonal. We have to make diagonal
392* and offdiagonal elements real, and store them in D and E.
393* For that, for real precision just copy the diag and offdiag
394* to D and E while for the COMPLEX case the bulge chasing is
395* performed to convert the hermetian tridiagonal to symmetric
396* tridiagonal. A simpler conversion formula might be used, but then
397* updating the Q matrix will be required and based if Q is generated
398* or not this might complicate the story.
399*
400 IF( kd.EQ.1 ) THEN
401 DO 50 i = 1, n
402 d( i ) = real( ab( abdpos, i ) )
403 50 CONTINUE
404*
405* make off-diagonal elements real and copy them to E
406*
407 IF( upper ) THEN
408 DO 60 i = 1, n - 1
409 tmp = ab( abofdpos, i+1 )
410 abstmp = abs( tmp )
411 ab( abofdpos, i+1 ) = abstmp
412 e( i ) = abstmp
413 IF( abstmp.NE.rzero ) THEN
414 tmp = tmp / abstmp
415 ELSE
416 tmp = one
417 END IF
418 IF( i.LT.n-1 )
419 $ ab( abofdpos, i+2 ) = ab( abofdpos, i+2 )*tmp
420C IF( WANTZ ) THEN
421C CALL CSCAL( N, CONJG( TMP ), Q( 1, I+1 ), 1 )
422C END IF
423 60 CONTINUE
424 ELSE
425 DO 70 i = 1, n - 1
426 tmp = ab( abofdpos, i )
427 abstmp = abs( tmp )
428 ab( abofdpos, i ) = abstmp
429 e( i ) = abstmp
430 IF( abstmp.NE.rzero ) THEN
431 tmp = tmp / abstmp
432 ELSE
433 tmp = one
434 END IF
435 IF( i.LT.n-1 )
436 $ ab( abofdpos, i+1 ) = ab( abofdpos, i+1 )*tmp
437C IF( WANTQ ) THEN
438C CALL CSCAL( N, TMP, Q( 1, I+1 ), 1 )
439C END IF
440 70 CONTINUE
441 ENDIF
442*
443 hous( 1 ) = 1
444 work( 1 ) = 1
445 RETURN
446 END IF
447*
448* Main code start here.
449* Reduce the hermitian band of A to a tridiagonal matrix.
450*
451 thgrsiz = n
452 grsiz = 1
453 shift = 3
454 nbtiles = ceiling( real(n)/real(kd) )
455 stepercol = ceiling( real(shift)/real(grsiz) )
456 thgrnb = ceiling( real(n-1)/real(thgrsiz) )
457*
458 CALL clacpy( "A", kd+1, n, ab, ldab, work( apos ), lda )
459 CALL claset( "A", kd, n, zero, zero, work( awpos ), lda )
460*
461*
462* openMP parallelisation start here
463*
464#if defined(_OPENMP)
465!$OMP PARALLEL PRIVATE( TID, THGRID, BLKLASTIND )
466!$OMP$ PRIVATE( THED, I, M, K, ST, ED, STT, SWEEPID )
467!$OMP$ PRIVATE( MYID, TTYPE, COLPT, STIND, EDIND )
468!$OMP$ SHARED ( UPLO, WANTQ, INDV, INDTAU, HOUS, WORK)
469!$OMP$ SHARED ( N, KD, IB, NBTILES, LDA, LDV, INDA )
470!$OMP$ SHARED ( STEPERCOL, THGRNB, THGRSIZ, GRSIZ, SHIFT )
471!$OMP MASTER
472#endif
473*
474* main bulge chasing loop
475*
476 DO 100 thgrid = 1, thgrnb
477 stt = (thgrid-1)*thgrsiz+1
478 thed = min( (stt + thgrsiz -1), (n-1))
479 DO 110 i = stt, n-1
480 ed = min( i, thed )
481 IF( stt.GT.ed ) EXIT
482 DO 120 m = 1, stepercol
483 st = stt
484 DO 130 sweepid = st, ed
485 DO 140 k = 1, grsiz
486 myid = (i-sweepid)*(stepercol*grsiz)
487 $ + (m-1)*grsiz + k
488 IF ( myid.EQ.1 ) THEN
489 ttype = 1
490 ELSE
491 ttype = mod( myid, 2 ) + 2
492 ENDIF
493
494 IF( ttype.EQ.2 ) THEN
495 colpt = (myid/2)*kd + sweepid
496 stind = colpt-kd+1
497 edind = min(colpt,n)
498 blklastind = colpt
499 ELSE
500 colpt = ((myid+1)/2)*kd + sweepid
501 stind = colpt-kd+1
502 edind = min(colpt,n)
503 IF( ( stind.GE.edind-1 ).AND.
504 $ ( edind.EQ.n ) ) THEN
505 blklastind = n
506 ELSE
507 blklastind = 0
508 ENDIF
509 ENDIF
510*
511* Call the kernel
512*
513#if defined(_OPENMP) && _OPENMP >= 201307
514 IF( ttype.NE.1 ) THEN
515!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
516!$OMP$ DEPEND(in:WORK(MYID-1))
517!$OMP$ DEPEND(out:WORK(MYID))
518 tid = omp_get_thread_num()
519 CALL chb2st_kernels( uplo, wantq, ttype,
520 $ stind, edind, sweepid, n, kd, ib,
521 $ work( inda ), lda,
522 $ hous( indv ), hous( indtau ), ldv,
523 $ work( indw + tid*kd ) )
524!$OMP END TASK
525 ELSE
526!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
527!$OMP$ DEPEND(out:WORK(MYID))
528 tid = omp_get_thread_num()
529 CALL chb2st_kernels( uplo, wantq, ttype,
530 $ stind, edind, sweepid, n, kd, ib,
531 $ work( inda ), lda,
532 $ hous( indv ), hous( indtau ), ldv,
533 $ work( indw + tid*kd ) )
534!$OMP END TASK
535 ENDIF
536#else
537 CALL chb2st_kernels( uplo, wantq, ttype,
538 $ stind, edind, sweepid, n, kd, ib,
539 $ work( inda ), lda,
540 $ hous( indv ), hous( indtau ), ldv,
541 $ work( indw ) )
542#endif
543 IF ( blklastind.GE.(n-1) ) THEN
544 stt = stt + 1
545 EXIT
546 ENDIF
547 140 CONTINUE
548 130 CONTINUE
549 120 CONTINUE
550 110 CONTINUE
551 100 CONTINUE
552*
553#if defined(_OPENMP)
554!$OMP END MASTER
555!$OMP END PARALLEL
556#endif
557*
558* Copy the diagonal from A to D. Note that D is REAL thus only
559* the Real part is needed, the imaginary part should be zero.
560*
561 DO 150 i = 1, n
562 d( i ) = real( work( dpos+(i-1)*lda ) )
563 150 CONTINUE
564*
565* Copy the off diagonal from A to E. Note that E is REAL thus only
566* the Real part is needed, the imaginary part should be zero.
567*
568 IF( upper ) THEN
569 DO 160 i = 1, n-1
570 e( i ) = real( work( ofdpos+i*lda ) )
571 160 CONTINUE
572 ELSE
573 DO 170 i = 1, n-1
574 e( i ) = real( work( ofdpos+(i-1)*lda ) )
575 170 CONTINUE
576 ENDIF
577*
578 hous( 1 ) = lhmin
579 work( 1 ) = sroundup_lwork(lwmin)
580 RETURN
581*
582* End of CHETRD_HB2ST
583*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine chb2st_kernels(uplo, wantz, ttype, st, ed, sweep, n, nb, ib, a, lda, v, tau, ldvt, work)
CHB2ST_KERNELS
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
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: