SRC/zhptrd.f

   1 *> \brief \b ZHPTRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHPTRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   D( * ), E( * )
  29 *       COMPLEX*16         AP( * ), TAU( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> ZHPTRD reduces a complex Hermitian matrix A stored in packed form to
  39 *> real symmetric tridiagonal form T by a unitary similarity
  40 *> transformation: Q**H * A * Q = T.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] UPLO
  47 *> \verbatim
  48 *>          UPLO is CHARACTER*1
  49 *>          = 'U':  Upper triangle of A is stored;
  50 *>          = 'L':  Lower triangle of A is stored.
  51 *> \endverbatim
  52 *>
  53 *> \param[in] N
  54 *> \verbatim
  55 *>          N is INTEGER
  56 *>          The order of the matrix A.  N >= 0.
  57 *> \endverbatim
  58 *>
  59 *> \param[in,out] AP
  60 *> \verbatim
  61 *>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
  62 *>          On entry, the upper or lower triangle of the Hermitian matrix
  63 *>          A, packed columnwise in a linear array.  The j-th column of A
  64 *>          is stored in the array AP as follows:
  65 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  66 *>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
  67 *>          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  68 *>          of A are overwritten by the corresponding elements of the
  69 *>          tridiagonal matrix T, and the elements above the first
  70 *>          superdiagonal, with the array TAU, represent the unitary
  71 *>          matrix Q as a product of elementary reflectors; if UPLO
  72 *>          = 'L', the diagonal and first subdiagonal of A are over-
  73 *>          written by the corresponding elements of the tridiagonal
  74 *>          matrix T, and the elements below the first subdiagonal, with
  75 *>          the array TAU, represent the unitary matrix Q as a product
  76 *>          of elementary reflectors. See Further Details.
  77 *> \endverbatim
  78 *>
  79 *> \param[out] D
  80 *> \verbatim
  81 *>          D is DOUBLE PRECISION array, dimension (N)
  82 *>          The diagonal elements of the tridiagonal matrix T:
  83 *>          D(i) = A(i,i).
  84 *> \endverbatim
  85 *>
  86 *> \param[out] E
  87 *> \verbatim
  88 *>          E is DOUBLE PRECISION array, dimension (N-1)
  89 *>          The off-diagonal elements of the tridiagonal matrix T:
  90 *>          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  91 *> \endverbatim
  92 *>
  93 *> \param[out] TAU
  94 *> \verbatim
  95 *>          TAU is COMPLEX*16 array, dimension (N-1)
  96 *>          The scalar factors of the elementary reflectors (see Further
  97 *>          Details).
  98 *> \endverbatim
  99 *>
 100 *> \param[out] INFO
 101 *> \verbatim
 102 *>          INFO is INTEGER
 103 *>          = 0:  successful exit
 104 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 105 *> \endverbatim
 106 *
 107 *  Authors:
 108 *  ========
 109 *
 110 *> \author Univ. of Tennessee
 111 *> \author Univ. of California Berkeley
 112 *> \author Univ. of Colorado Denver
 113 *> \author NAG Ltd.
 114 *
 115 *> \date November 2011
 116 *
 117 *> \ingroup complex16OTHERcomputational
 118 *
 119 *> \par Further Details:
 120 *  =====================
 121 *>
 122 *> \verbatim
 123 *>
 124 *>  If UPLO = 'U', the matrix Q is represented as a product of elementary
 125 *>  reflectors
 126 *>
 127 *>     Q = H(n-1) . . . H(2) H(1).
 128 *>
 129 *>  Each H(i) has the form
 130 *>
 131 *>     H(i) = I - tau * v * v**H
 132 *>
 133 *>  where tau is a complex scalar, and v is a complex vector with
 134 *>  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
 135 *>  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
 136 *>
 137 *>  If UPLO = 'L', the matrix Q is represented as a product of elementary
 138 *>  reflectors
 139 *>
 140 *>     Q = H(1) H(2) . . . H(n-1).
 141 *>
 142 *>  Each H(i) has the form
 143 *>
 144 *>     H(i) = I - tau * v * v**H
 145 *>
 146 *>  where tau is a complex scalar, and v is a complex vector with
 147 *>  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
 148 *>  overwriting A(i+2:n,i), and tau is stored in TAU(i).
 149 *> \endverbatim
 150 *>
 151 *  =====================================================================
 152       SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
 153 *
 154 *  -- LAPACK computational routine (version 3.4.0) --
 155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 157 *     November 2011
 158 *
 159 *     .. Scalar Arguments ..
 160       CHARACTER          UPLO
 161       INTEGER            INFO, N
 162 *     ..
 163 *     .. Array Arguments ..
 164       DOUBLE PRECISION   D( * ), E( * )
 165       COMPLEX*16         AP( * ), TAU( * )
 166 *     ..
 167 *
 168 *  =====================================================================
 169 *
 170 *     .. Parameters ..
 171       COMPLEX*16         ONE, ZERO, HALF
 172       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
 173      $                   ZERO = ( 0.0D+0, 0.0D+0 ),
 174      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
 175 *     ..
 176 *     .. Local Scalars ..
 177       LOGICAL            UPPER
 178       INTEGER            I, I1, I1I1, II
 179       COMPLEX*16         ALPHA, TAUI
 180 *     ..
 181 *     .. External Subroutines ..
 182       EXTERNAL           XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG
 183 *     ..
 184 *     .. External Functions ..
 185       LOGICAL            LSAME
 186       COMPLEX*16         ZDOTC
 187       EXTERNAL           LSAME, ZDOTC
 188 *     ..
 189 *     .. Intrinsic Functions ..
 190       INTRINSIC          DBLE
 191 *     ..
 192 *     .. Executable Statements ..
 193 *
 194 *     Test the input parameters
 195 *
 196       INFO = 0
 197       UPPER = LSAME( UPLO, 'U' )
 198       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 199          INFO = -1
 200       ELSE IF( N.LT.0 ) THEN
 201          INFO = -2
 202       END IF
 203       IF( INFO.NE.0 ) THEN
 204          CALL XERBLA( 'ZHPTRD', -INFO )
 205          RETURN
 206       END IF
 207 *
 208 *     Quick return if possible
 209 *
 210       IF( N.LE.0 )
 211      $   RETURN
 212 *
 213       IF( UPPER ) THEN
 214 *
 215 *        Reduce the upper triangle of A.
 216 *        I1 is the index in AP of A(1,I+1).
 217 *
 218          I1 = N*( N-1 ) / 2 + 1
 219          AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) )
 220          DO 10 I = N - 1, 1, -1
 221 *
 222 *           Generate elementary reflector H(i) = I - tau * v * v**H
 223 *           to annihilate A(1:i-1,i+1)
 224 *
 225             ALPHA = AP( I1+I-1 )
 226             CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI )
 227             E( I ) = ALPHA
 228 *
 229             IF( TAUI.NE.ZERO ) THEN
 230 *
 231 *              Apply H(i) from both sides to A(1:i,1:i)
 232 *
 233                AP( I1+I-1 ) = ONE
 234 *
 235 *              Compute  y := tau * A * v  storing y in TAU(1:i)
 236 *
 237                CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
 238      $                     1 )
 239 *
 240 *              Compute  w := y - 1/2 * tau * (y**H *v) * v
 241 *
 242                ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 )
 243                CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
 244 *
 245 *              Apply the transformation as a rank-2 update:
 246 *                 A := A - v * w**H - w * v**H
 247 *
 248                CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
 249 *
 250             END IF
 251             AP( I1+I-1 ) = E( I )
 252             D( I+1 ) = AP( I1+I )
 253             TAU( I ) = TAUI
 254             I1 = I1 - I
 255    10    CONTINUE
 256          D( 1 ) = AP( 1 )
 257       ELSE
 258 *
 259 *        Reduce the lower triangle of A. II is the index in AP of
 260 *        A(i,i) and I1I1 is the index of A(i+1,i+1).
 261 *
 262          II = 1
 263          AP( 1 ) = DBLE( AP( 1 ) )
 264          DO 20 I = 1, N - 1
 265             I1I1 = II + N - I + 1
 266 *
 267 *           Generate elementary reflector H(i) = I - tau * v * v**H
 268 *           to annihilate A(i+2:n,i)
 269 *
 270             ALPHA = AP( II+1 )
 271             CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI )
 272             E( I ) = ALPHA
 273 *
 274             IF( TAUI.NE.ZERO ) THEN
 275 *
 276 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
 277 *
 278                AP( II+1 ) = ONE
 279 *
 280 *              Compute  y := tau * A * v  storing y in TAU(i:n-1)
 281 *
 282                CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
 283      $                     ZERO, TAU( I ), 1 )
 284 *
 285 *              Compute  w := y - 1/2 * tau * (y**H *v) * v
 286 *
 287                ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ),
 288      $                 1 )
 289                CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
 290 *
 291 *              Apply the transformation as a rank-2 update:
 292 *                 A := A - v * w**H - w * v**H
 293 *
 294                CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
 295      $                     AP( I1I1 ) )
 296 *
 297             END IF
 298             AP( II+1 ) = E( I )
 299             D( I ) = AP( II )
 300             TAU( I ) = TAUI
 301             II = I1I1
 302    20    CONTINUE
 303          D( N ) = AP( II )
 304       END IF
 305 *
 306       RETURN
 307 *
 308 *     End of ZHPTRD
 309 *
 310       END