chbtrd.f

Go to the documentation of this file.

*> \brief \b CHBTRD
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CHBTRD + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbtrd.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbtrd.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbtrd.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
*                          WORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          UPLO, VECT
*       INTEGER            INFO, KD, LDAB, LDQ, N
*       ..
*       .. Array Arguments ..
*       REAL               D( * ), E( * )
*       COMPLEX            AB( LDAB, * ), Q( LDQ, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CHBTRD reduces a complex Hermitian band matrix A to real symmetric
*> tridiagonal form T by a unitary similarity transformation:
*> Q**H * A * Q = T.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] VECT
*> \verbatim
*>          VECT is CHARACTER*1
*>          = 'N':  do not form Q;
*>          = 'V':  form Q;
*>          = 'U':  update a matrix X, by forming X*Q.
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  Upper triangle of A is stored;
*>          = 'L':  Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in] KD
*> \verbatim
*>          KD is INTEGER
*>          The number of superdiagonals of the matrix A if UPLO = 'U',
*>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
*> \endverbatim
*>
*> \param[in,out] AB
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDAB
*> \verbatim
*>          LDAB is INTEGER
*>          The leading dimension of the array AB.  LDAB >= KD+1.
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*>          D is REAL array, dimension (N)
*>          The diagonal elements of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[out] E
*> \verbatim
*>          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'.
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          Q is COMPLEX array, dimension (LDQ,N)
*>          On entry, if VECT = 'U', then Q must contain an N-by-N
*>          matrix X; if VECT = 'N' or 'V', then Q need not be set.
*>
*>          On exit:
*>          if VECT = 'V', Q contains the N-by-N unitary matrix Q;
*>          if VECT = 'U', Q contains the product X*Q;
*>          if VECT = 'N', the array Q is not referenced.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.
*>          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX array, dimension (N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hbtrd
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  Modified by Linda Kaufman, Bell Labs.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE chbtrd( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
     $                   WORK, INFO )
*
*  -- 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, VECT
      INTEGER            INFO, KD, LDAB, LDQ, N
*     ..
*     .. Array Arguments ..
      REAL               D( * ), E( * )
      COMPLEX            AB( LDAB, * ), Q( LDQ, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      parameter( zero = 0.0e+0 )
      COMPLEX            CZERO, CONE
      parameter( czero = ( 0.0e+0, 0.0e+0 ),
     $                   cone = ( 1.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            INITQ, UPPER, WANTQ
      INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
     $                   j1, j1end, j1inc, j2, jend, jin, jinc, k, kd1,
     $                   kdm1, kdn, l, last, lend, nq, nr, nrt
      REAL               ABST
      COMPLEX            T, TEMP
*     ..
*     .. External Subroutines ..
      EXTERNAL           clacgv, clar2v, clargv, clartg, clartv, claset,
     $                   crot, cscal, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, conjg, max, min, real
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      initq = lsame( vect, 'V' )
      wantq = initq .OR. lsame( vect, 'U' )
      upper = lsame( uplo, 'U' )
      kd1 = kd + 1
      kdm1 = kd - 1
      incx = ldab - 1
      iqend = 1
*
      info = 0
      IF( .NOT.wantq .AND. .NOT.lsame( vect, 'N' ) ) THEN
         info = -1
      ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( kd.LT.0 ) THEN
         info = -4
      ELSE IF( ldab.LT.kd1 ) THEN
         info = -6
      ELSE IF( ldq.LT.max( 1, n ) .AND. wantq ) THEN
         info = -10
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CHBTRD', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
*
*     Initialize Q to the unit matrix, if needed
*
      IF( initq )
     $   CALL claset( 'Full', n, n, czero, cone, q, ldq )
*
*     Wherever possible, plane rotations are generated and applied in
*     vector operations of length NR over the index set J1:J2:KD1.
*
*     The real cosines and complex sines of the plane rotations are
*     stored in the arrays D and WORK.
*
      inca = kd1*ldab
      kdn = min( n-1, kd )
      IF( upper ) THEN
*
         IF( kd.GT.1 ) THEN
*
*           Reduce to complex Hermitian tridiagonal form, working with
*           the upper triangle
*
            nr = 0
            j1 = kdn + 2
            j2 = 1
*
            ab( kd1, 1 ) = real( ab( kd1, 1 ) )
            DO 90 i = 1, n - 2
*
*              Reduce i-th row of matrix to tridiagonal form
*
               DO 80 k = kdn + 1, 2, -1
                  j1 = j1 + kdn
                  j2 = j2 + kdn
*
                  IF( nr.GT.0 ) THEN
*
*                    generate plane rotations to annihilate nonzero
*                    elements which have been created outside the band
*
                     CALL clargv( nr, ab( 1, j1-1 ), inca, work( j1 ),
     $                            kd1, d( j1 ), kd1 )
*
*                    apply rotations from the right
*
*
*                    Dependent on the the number of diagonals either
*                    CLARTV or CROT is used
*
                     IF( nr.GE.2*kd-1 ) THEN
                        DO 10 l = 1, kd - 1
                           CALL clartv( nr, ab( l+1, j1-1 ), inca,
     $                                  ab( l, j1 ), inca, d( j1 ),
     $                                  work( j1 ), kd1 )
   10                   CONTINUE
*
                     ELSE
                        jend = j1 + ( nr-1 )*kd1
                        DO 20 jinc = j1, jend, kd1
                           CALL crot( kdm1, ab( 2, jinc-1 ), 1,
     $                                ab( 1, jinc ), 1, d( jinc ),
     $                                work( jinc ) )
   20                   CONTINUE
                     END IF
                  END IF
*
*
                  IF( k.GT.2 ) THEN
                     IF( k.LE.n-i+1 ) THEN
*
*                       generate plane rotation to annihilate a(i,i+k-1)
*                       within the band
*
                        CALL clartg( ab( kd-k+3, i+k-2 ),
     $                               ab( kd-k+2, i+k-1 ), d( i+k-1 ),
     $                               work( i+k-1 ), temp )
                        ab( kd-k+3, i+k-2 ) = temp
*
*                       apply rotation from the right
*
                        CALL crot( k-3, ab( kd-k+4, i+k-2 ), 1,
     $                             ab( kd-k+3, i+k-1 ), 1, d( i+k-1 ),
     $                             work( i+k-1 ) )
                     END IF
                     nr = nr + 1
                     j1 = j1 - kdn - 1
                  END IF
*
*                 apply plane rotations from both sides to diagonal
*                 blocks
*
                  IF( nr.GT.0 )
     $               CALL clar2v( nr, ab( kd1, j1-1 ), ab( kd1, j1 ),
     $                            ab( kd, j1 ), inca, d( j1 ),
     $                            work( j1 ), kd1 )
*
*                 apply plane rotations from the left
*
                  IF( nr.GT.0 ) THEN
                     CALL clacgv( nr, work( j1 ), kd1 )
                     IF( 2*kd-1.LT.nr ) THEN
*
*                    Dependent on the the number of diagonals either
*                    CLARTV or CROT is used
*
                        DO 30 l = 1, kd - 1
                           IF( j2+l.GT.n ) THEN
                              nrt = nr - 1
                           ELSE
                              nrt = nr
                           END IF
                           IF( nrt.GT.0 )
     $                        CALL clartv( nrt, ab( kd-l, j1+l ), inca,
     $                                     ab( kd-l+1, j1+l ), inca,
     $                                     d( j1 ), work( j1 ), kd1 )
   30                   CONTINUE
                     ELSE
                        j1end = j1 + kd1*( nr-2 )
                        IF( j1end.GE.j1 ) THEN
                           DO 40 jin = j1, j1end, kd1
                              CALL crot( kd-1, ab( kd-1, jin+1 ), incx,
     $                                   ab( kd, jin+1 ), incx,
     $                                   d( jin ), work( jin ) )
   40                      CONTINUE
                        END IF
                        lend = min( kdm1, n-j2 )
                        last = j1end + kd1
                        IF( lend.GT.0 )
     $                     CALL crot( lend, ab( kd-1, last+1 ), incx,
     $                                ab( kd, last+1 ), incx, d( last ),
     $                                work( last ) )
                     END IF
                  END IF
*
                  IF( wantq ) THEN
*
*                    accumulate product of plane rotations in Q
*
                     IF( initq ) THEN
*
*                 take advantage of the fact that Q was
*                 initially the Identity matrix
*
                        iqend = max( iqend, j2 )
                        i2 = max( 0, k-3 )
                        iqaend = 1 + i*kd
                        IF( k.EQ.2 )
     $                     iqaend = iqaend + kd
                        iqaend = min( iqaend, iqend )
                        DO 50 j = j1, j2, kd1
                           ibl = i - i2 / kdm1
                           i2 = i2 + 1
                           iqb = max( 1, j-ibl )
                           nq = 1 + iqaend - iqb
                           iqaend = min( iqaend+kd, iqend )
                           CALL crot( nq, q( iqb, j-1 ), 1, q( iqb, j ),
     $                                1, d( j ), conjg( work( j ) ) )
   50                   CONTINUE
                     ELSE
*
                        DO 60 j = j1, j2, kd1
                           CALL crot( n, q( 1, j-1 ), 1, q( 1, j ), 1,
     $                                d( j ), conjg( work( j ) ) )
   60                   CONTINUE
                     END IF
*
                  END IF
*
                  IF( j2+kdn.GT.n ) THEN
*
*                    adjust J2 to keep within the bounds of the matrix
*
                     nr = nr - 1
                     j2 = j2 - kdn - 1
                  END IF
*
                  DO 70 j = j1, j2, kd1
*
*                    create nonzero element a(j-1,j+kd) outside the band
*                    and store it in WORK
*
                     work( j+kd ) = work( j )*ab( 1, j+kd )
                     ab( 1, j+kd ) = d( j )*ab( 1, j+kd )
   70             CONTINUE
   80          CONTINUE
   90       CONTINUE
         END IF
*
         IF( kd.GT.0 ) THEN
*
*           make off-diagonal elements real and copy them to E
*
            DO 100 i = 1, n - 1
               t = ab( kd, i+1 )
               abst = abs( t )
               ab( kd, i+1 ) = abst
               e( i ) = abst
               IF( abst.NE.zero ) THEN
                  t = t / abst
               ELSE
                  t = cone
               END IF
               IF( i.LT.n-1 )
     $            ab( kd, i+2 ) = ab( kd, i+2 )*t
               IF( wantq ) THEN
                  CALL cscal( n, conjg( t ), q( 1, i+1 ), 1 )
               END IF
  100       CONTINUE
         ELSE
*
*           set E to zero if original matrix was diagonal
*
            DO 110 i = 1, n - 1
               e( i ) = zero
  110       CONTINUE
         END IF
*
*        copy diagonal elements to D
*
         DO 120 i = 1, n
            d( i ) = real( ab( kd1, i ) )
  120    CONTINUE
*
      ELSE
*
         IF( kd.GT.1 ) THEN
*
*           Reduce to complex Hermitian tridiagonal form, working with
*           the lower triangle
*
            nr = 0
            j1 = kdn + 2
            j2 = 1
*
            ab( 1, 1 ) = real( ab( 1, 1 ) )
            DO 210 i = 1, n - 2
*
*              Reduce i-th column of matrix to tridiagonal form
*
               DO 200 k = kdn + 1, 2, -1
                  j1 = j1 + kdn
                  j2 = j2 + kdn
*
                  IF( nr.GT.0 ) THEN
*
*                    generate plane rotations to annihilate nonzero
*                    elements which have been created outside the band
*
                     CALL clargv( nr, ab( kd1, j1-kd1 ), inca,
     $                            work( j1 ), kd1, d( j1 ), kd1 )
*
*                    apply plane rotations from one side
*
*
*                    Dependent on the the number of diagonals either
*                    CLARTV or CROT is used
*
                     IF( nr.GT.2*kd-1 ) THEN
                        DO 130 l = 1, kd - 1
                           CALL clartv( nr, ab( kd1-l, j1-kd1+l ), inca,
     $                                  ab( kd1-l+1, j1-kd1+l ), inca,
     $                                  d( j1 ), work( j1 ), kd1 )
  130                   CONTINUE
                     ELSE
                        jend = j1 + kd1*( nr-1 )
                        DO 140 jinc = j1, jend, kd1
                           CALL crot( kdm1, ab( kd, jinc-kd ), incx,
     $                                ab( kd1, jinc-kd ), incx,
     $                                d( jinc ), work( jinc ) )
  140                   CONTINUE
                     END IF
*
                  END IF
*
                  IF( k.GT.2 ) THEN
                     IF( k.LE.n-i+1 ) THEN
*
*                       generate plane rotation to annihilate a(i+k-1,i)
*                       within the band
*
                        CALL clartg( ab( k-1, i ), ab( k, i ),
     $                               d( i+k-1 ), work( i+k-1 ), temp )
                        ab( k-1, i ) = temp
*
*                       apply rotation from the left
*
                        CALL crot( k-3, ab( k-2, i+1 ), ldab-1,
     $                             ab( k-1, i+1 ), ldab-1, d( i+k-1 ),
     $                             work( i+k-1 ) )
                     END IF
                     nr = nr + 1
                     j1 = j1 - kdn - 1
                  END IF
*
*                 apply plane rotations from both sides to diagonal
*                 blocks
*
                  IF( nr.GT.0 )
     $               CALL clar2v( nr, ab( 1, j1-1 ), ab( 1, j1 ),
     $                            ab( 2, j1-1 ), inca, d( j1 ),
     $                            work( j1 ), kd1 )
*
*                 apply plane rotations from the right
*
*
*                    Dependent on the the number of diagonals either
*                    CLARTV or CROT is used
*
                  IF( nr.GT.0 ) THEN
                     CALL clacgv( nr, work( j1 ), kd1 )
                     IF( nr.GT.2*kd-1 ) THEN
                        DO 150 l = 1, kd - 1
                           IF( j2+l.GT.n ) THEN
                              nrt = nr - 1
                           ELSE
                              nrt = nr
                           END IF
                           IF( nrt.GT.0 )
     $                        CALL clartv( nrt, ab( l+2, j1-1 ), inca,
     $                                     ab( l+1, j1 ), inca, d( j1 ),
     $                                     work( j1 ), kd1 )
  150                   CONTINUE
                     ELSE
                        j1end = j1 + kd1*( nr-2 )
                        IF( j1end.GE.j1 ) THEN
                           DO 160 j1inc = j1, j1end, kd1
                              CALL crot( kdm1, ab( 3, j1inc-1 ), 1,
     $                                   ab( 2, j1inc ), 1, d( j1inc ),
     $                                   work( j1inc ) )
  160                      CONTINUE
                        END IF
                        lend = min( kdm1, n-j2 )
                        last = j1end + kd1
                        IF( lend.GT.0 )
     $                     CALL crot( lend, ab( 3, last-1 ), 1,
     $                                ab( 2, last ), 1, d( last ),
     $                                work( last ) )
                     END IF
                  END IF
*
*
*
                  IF( wantq ) THEN
*
*                    accumulate product of plane rotations in Q
*
                     IF( initq ) THEN
*
*                 take advantage of the fact that Q was
*                 initially the Identity matrix
*
                        iqend = max( iqend, j2 )
                        i2 = max( 0, k-3 )
                        iqaend = 1 + i*kd
                        IF( k.EQ.2 )
     $                     iqaend = iqaend + kd
                        iqaend = min( iqaend, iqend )
                        DO 170 j = j1, j2, kd1
                           ibl = i - i2 / kdm1
                           i2 = i2 + 1
                           iqb = max( 1, j-ibl )
                           nq = 1 + iqaend - iqb
                           iqaend = min( iqaend+kd, iqend )
                           CALL crot( nq, q( iqb, j-1 ), 1, q( iqb, j ),
     $                                1, d( j ), work( j ) )
  170                   CONTINUE
                     ELSE
*
                        DO 180 j = j1, j2, kd1
                           CALL crot( n, q( 1, j-1 ), 1, q( 1, j ), 1,
     $                                d( j ), work( j ) )
  180                   CONTINUE
                     END IF
                  END IF
*
                  IF( j2+kdn.GT.n ) THEN
*
*                    adjust J2 to keep within the bounds of the matrix
*
                     nr = nr - 1
                     j2 = j2 - kdn - 1
                  END IF
*
                  DO 190 j = j1, j2, kd1
*
*                    create nonzero element a(j+kd,j-1) outside the
*                    band and store it in WORK
*
                     work( j+kd ) = work( j )*ab( kd1, j )
                     ab( kd1, j ) = d( j )*ab( kd1, j )
  190             CONTINUE
  200          CONTINUE
  210       CONTINUE
         END IF
*
         IF( kd.GT.0 ) THEN
*
*           make off-diagonal elements real and copy them to E
*
            DO 220 i = 1, n - 1
               t = ab( 2, i )
               abst = abs( t )
               ab( 2, i ) = abst
               e( i ) = abst
               IF( abst.NE.zero ) THEN
                  t = t / abst
               ELSE
                  t = cone
               END IF
               IF( i.LT.n-1 )
     $            ab( 2, i+1 ) = ab( 2, i+1 )*t
               IF( wantq ) THEN
                  CALL cscal( n, t, q( 1, i+1 ), 1 )
               END IF
  220       CONTINUE
         ELSE
*
*           set E to zero if original matrix was diagonal
*
            DO 230 i = 1, n - 1
               e( i ) = zero
  230       CONTINUE
         END IF
*
*        copy diagonal elements to D
*
         DO 240 i = 1, n
            d( i ) = real( ab( 1, i ) )
  240    CONTINUE
      END IF
*
      RETURN
*
*     End of CHBTRD
*
      END