zhb2st_kernels.f

Go to the documentation of this file.

*> \brief \b ZHB2ST_KERNELS
*
*  @precisions fortran z -> s d c
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download ZHB2ST_KERNELS + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhb2st_kernels.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhb2st_kernels.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhb2st_kernels.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*       SUBROUTINE  ZHB2ST_KERNELS( UPLO, WANTZ, TTYPE,
*                                   ST, ED, SWEEP, N, NB, IB,
*                                   A, LDA, V, TAU, LDVT, WORK)
*
*       IMPLICIT NONE
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO
*       LOGICAL            WANTZ
*       INTEGER            TTYPE, ST, ED, SWEEP, N, NB, IB, LDA, LDVT
*       ..
*       .. Array Arguments ..
*       COMPLEX*16         A( LDA, * ), V( * ),
*                          TAU( * ), WORK( * )
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZHB2ST_KERNELS is an internal routine used by the ZHETRD_HB2ST
*> subroutine.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*> \endverbatim
*>
*> \param[in] WANTZ
*> \verbatim
*>          WANTZ is LOGICAL which indicate if Eigenvalue are requested or both
*>          Eigenvalue/Eigenvectors.
*> \endverbatim
*>
*> \param[in] TTYPE
*> \verbatim
*>          TTYPE is INTEGER
*> \endverbatim
*>
*> \param[in] ST
*> \verbatim
*>          ST is INTEGER
*>          internal parameter for indices.
*> \endverbatim
*>
*> \param[in] ED
*> \verbatim
*>          ED is INTEGER
*>          internal parameter for indices.
*> \endverbatim
*>
*> \param[in] SWEEP
*> \verbatim
*>          SWEEP is INTEGER
*>          internal parameter for indices.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER. The order of the matrix A.
*> \endverbatim
*>
*> \param[in] NB
*> \verbatim
*>          NB is INTEGER. The size of the band.
*> \endverbatim
*>
*> \param[in] IB
*> \verbatim
*>          IB is INTEGER.
*> \endverbatim
*>
*> \param[in, out] A
*> \verbatim
*>          A is COMPLEX*16 array. A pointer to the matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER. The leading dimension of the matrix A.
*> \endverbatim
*>
*> \param[out] V
*> \verbatim
*>          V is COMPLEX*16 array, dimension 2*n if eigenvalues only are
*>          requested or to be queried for vectors.
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*>          TAU is COMPLEX*16 array, dimension (2*n).
*>          The scalar factors of the Householder reflectors are stored
*>          in this array.
*> \endverbatim
*>
*> \param[in] LDVT
*> \verbatim
*>          LDVT is INTEGER.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array. Workspace of size nb.
*> \endverbatim
*>
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  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
*>
*> \endverbatim
*>
*> \ingroup hb2st_kernels
*>
*  =====================================================================
      SUBROUTINE  zhb2st_kernels( UPLO, WANTZ, TTYPE,
     $                            ST, ED, SWEEP, N, NB, IB,
     $                            A, LDA, V, TAU, LDVT, WORK)
*
      IMPLICIT NONE
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      LOGICAL            WANTZ
      INTEGER            TTYPE, ST, ED, SWEEP, N, NB, IB, LDA, LDVT
*     ..
*     .. Array Arguments ..
      COMPLEX*16         A( LDA, * ), V( * ),
     $                   TAU( * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         ZERO, ONE
      PARAMETER          ( ZERO = ( 0.0d+0, 0.0d+0 ),
     $                   one = ( 1.0d+0, 0.0d+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      INTEGER            I, J1, J2, LM, LN, VPOS, TAUPOS,
     $                   dpos, ofdpos, ajeter
      COMPLEX*16         CTMP
*     ..
*     .. External Subroutines ..
      EXTERNAL           zlarfg, zlarfx, zlarfy
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dconjg, mod
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     ..
*     .. Executable Statements ..
*
      ajeter = ib + ldvt
      upper = lsame( uplo, 'U' )
 
      IF( upper ) THEN
          dpos    = 2 * nb + 1
          ofdpos  = 2 * nb
      ELSE
          dpos    = 1
          ofdpos  = 2
      ENDIF
 
*
*     Upper case
*
      IF( upper ) THEN
*
          IF( wantz ) THEN
              vpos   = mod( sweep-1, 2 ) * n + st
              taupos = mod( sweep-1, 2 ) * n + st
          ELSE
              vpos   = mod( sweep-1, 2 ) * n + st
              taupos = mod( sweep-1, 2 ) * n + st
          ENDIF
*
          IF( ttype.EQ.1 ) THEN
              lm = ed - st + 1
*
              v( vpos ) = one
              DO 10 i = 1, lm-1
                  v( vpos+i )         = dconjg( a( ofdpos-i, st+i ) )
                  a( ofdpos-i, st+i ) = zero
   10         CONTINUE
              ctmp = dconjg( a( ofdpos, st ) )
              CALL zlarfg( lm, ctmp, v( vpos+1 ), 1,
     $                                       tau( taupos ) )
              a( ofdpos, st ) = ctmp
*
              lm = ed - st + 1
              CALL zlarfy( uplo, lm, v( vpos ), 1,
     $                     dconjg( tau( taupos ) ),
     $                     a( dpos, st ), lda-1, work)
          ENDIF
*
          IF( ttype.EQ.3 ) THEN
*
              lm = ed - st + 1
              CALL zlarfy( uplo, lm, v( vpos ), 1,
     $                     dconjg( tau( taupos ) ),
     $                     a( dpos, st ), lda-1, work)
          ENDIF
*
          IF( ttype.EQ.2 ) THEN
              j1 = ed+1
              j2 = min( ed+nb, n )
              ln = ed-st+1
              lm = j2-j1+1
              IF( lm.GT.0) THEN
                  CALL zlarfx( 'Left', ln, lm, v( vpos ),
     $                         dconjg( tau( taupos ) ),
     $                         a( dpos-nb, j1 ), lda-1, work)
*
                  IF( wantz ) THEN
                      vpos   = mod( sweep-1, 2 ) * n + j1
                      taupos = mod( sweep-1, 2 ) * n + j1
                  ELSE
                      vpos   = mod( sweep-1, 2 ) * n + j1
                      taupos = mod( sweep-1, 2 ) * n + j1
                  ENDIF
*
                  v( vpos ) = one
                  DO 30 i = 1, lm-1
                      v( vpos+i )          =
     $                                    dconjg( a( dpos-nb-i, j1+i ) )
                      a( dpos-nb-i, j1+i ) = zero
   30             CONTINUE
                  ctmp = dconjg( a( dpos-nb, j1 ) )
                  CALL zlarfg( lm, ctmp, v( vpos+1 ), 1,
     $                         tau( taupos ) )
                  a( dpos-nb, j1 ) = ctmp
*
                  CALL zlarfx( 'Right', ln-1, lm, v( vpos ),
     $                         tau( taupos ),
     $                         a( dpos-nb+1, j1 ), lda-1, work)
              ENDIF
          ENDIF
*
*     Lower case
*
      ELSE
*
          IF( wantz ) THEN
              vpos   = mod( sweep-1, 2 ) * n + st
              taupos = mod( sweep-1, 2 ) * n + st
          ELSE
              vpos   = mod( sweep-1, 2 ) * n + st
              taupos = mod( sweep-1, 2 ) * n + st
          ENDIF
*
          IF( ttype.EQ.1 ) THEN
              lm = ed - st + 1
*
              v( vpos ) = one
              DO 20 i = 1, lm-1
                  v( vpos+i )         = a( ofdpos+i, st-1 )
                  a( ofdpos+i, st-1 ) = zero
   20         CONTINUE
              CALL zlarfg( lm, a( ofdpos, st-1 ), v( vpos+1 ), 1,
     $                                       tau( taupos ) )
*
              lm = ed - st + 1
*
              CALL zlarfy( uplo, lm, v( vpos ), 1,
     $                     dconjg( tau( taupos ) ),
     $                     a( dpos, st ), lda-1, work)
 
          ENDIF
*
          IF( ttype.EQ.3 ) THEN
              lm = ed - st + 1
*
              CALL zlarfy( uplo, lm, v( vpos ), 1,
     $                     dconjg( tau( taupos ) ),
     $                     a( dpos, st ), lda-1, work)
 
          ENDIF
*
          IF( ttype.EQ.2 ) THEN
              j1 = ed+1
              j2 = min( ed+nb, n )
              ln = ed-st+1
              lm = j2-j1+1
*
              IF( lm.GT.0) THEN
                  CALL zlarfx( 'Right', lm, ln, v( vpos ),
     $                         tau( taupos ), a( dpos+nb, st ),
     $                         lda-1, work)
*
                  IF( wantz ) THEN
                      vpos   = mod( sweep-1, 2 ) * n + j1
                      taupos = mod( sweep-1, 2 ) * n + j1
                  ELSE
                      vpos   = mod( sweep-1, 2 ) * n + j1
                      taupos = mod( sweep-1, 2 ) * n + j1
                  ENDIF
*
                  v( vpos ) = one
                  DO 40 i = 1, lm-1
                      v( vpos+i )        = a( dpos+nb+i, st )
                      a( dpos+nb+i, st ) = zero
   40             CONTINUE
                  CALL zlarfg( lm, a( dpos+nb, st ), v( vpos+1 ), 1,
     $                                        tau( taupos ) )
*
                  CALL zlarfx( 'Left', lm, ln-1, v( vpos ),
     $                         dconjg( tau( taupos ) ),
     $                         a( dpos+nb-1, st+1 ), lda-1, work)
 
              ENDIF
          ENDIF
      ENDIF
*
      RETURN
*
*     End of ZHB2ST_KERNELS
*
      SUBROUTINE  zhb2st_kernels( UPLO, WANTZ, TTYPE, …
      END