src/LAPACK/zhptrd.f

   1       SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          UPLO
   9       INTEGER            INFO, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   D( * ), E( * )
  13       COMPLEX*16         AP( * ), TAU( * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  ZHPTRD reduces a complex Hermitian matrix A stored in packed form to
  20 *  real symmetric tridiagonal form T by a unitary similarity
  21 *  transformation: Q**H * A * Q = T.
  22 *
  23 *  Arguments
  24 *  =========
  25 *
  26 *  UPLO    (input) CHARACTER*1
  27 *          = 'U':  Upper triangle of A is stored;
  28 *          = 'L':  Lower triangle of A is stored.
  29 *
  30 *  N       (input) INTEGER
  31 *          The order of the matrix A.  N >= 0.
  32 *
  33 *  AP      (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
  34 *          On entry, the upper or lower triangle of the Hermitian matrix
  35 *          A, packed columnwise in a linear array.  The j-th column of A
  36 *          is stored in the array AP as follows:
  37 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  38 *          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
  39 *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  40 *          of A are overwritten by the corresponding elements of the
  41 *          tridiagonal matrix T, and the elements above the first
  42 *          superdiagonal, with the array TAU, represent the unitary
  43 *          matrix Q as a product of elementary reflectors; if UPLO
  44 *          = 'L', the diagonal and first subdiagonal of A are over-
  45 *          written by the corresponding elements of the tridiagonal
  46 *          matrix T, and the elements below the first subdiagonal, with
  47 *          the array TAU, represent the unitary matrix Q as a product
  48 *          of elementary reflectors. See Further Details.
  49 *
  50 *  D       (output) DOUBLE PRECISION array, dimension (N)
  51 *          The diagonal elements of the tridiagonal matrix T:
  52 *          D(i) = A(i,i).
  53 *
  54 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
  55 *          The off-diagonal elements of the tridiagonal matrix T:
  56 *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  57 *
  58 *  TAU     (output) COMPLEX*16 array, dimension (N-1)
  59 *          The scalar factors of the elementary reflectors (see Further
  60 *          Details).
  61 *
  62 *  INFO    (output) INTEGER
  63 *          = 0:  successful exit
  64 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  65 *
  66 *  Further Details
  67 *  ===============
  68 *
  69 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
  70 *  reflectors
  71 *
  72 *     Q = H(n-1) . . . H(2) H(1).
  73 *
  74 *  Each H(i) has the form
  75 *
  76 *     H(i) = I - tau * v * v'
  77 *
  78 *  where tau is a complex scalar, and v is a complex vector with
  79 *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
  80 *  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
  81 *
  82 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
  83 *  reflectors
  84 *
  85 *     Q = H(1) H(2) . . . H(n-1).
  86 *
  87 *  Each H(i) has the form
  88 *
  89 *     H(i) = I - tau * v * v'
  90 *
  91 *  where tau is a complex scalar, and v is a complex vector with
  92 *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
  93 *  overwriting A(i+2:n,i), and tau is stored in TAU(i).
  94 *
  95 *  =====================================================================
  96 *
  97 *     .. Parameters ..
  98       COMPLEX*16         ONE, ZERO, HALF
  99       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
 100      $                   ZERO = ( 0.0D+0, 0.0D+0 ),
 101      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
 102 *     ..
 103 *     .. Local Scalars ..
 104       LOGICAL            UPPER
 105       INTEGER            I, I1, I1I1, II
 106       COMPLEX*16         ALPHA, TAUI
 107 *     ..
 108 *     .. External Subroutines ..
 109       EXTERNAL           XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG
 110 *     ..
 111 *     .. External Functions ..
 112       LOGICAL            LSAME
 113       COMPLEX*16         ZDOTC
 114       EXTERNAL           LSAME, ZDOTC
 115 *     ..
 116 *     .. Intrinsic Functions ..
 117       INTRINSIC          DBLE
 118 *     ..
 119 *     .. Executable Statements ..
 120 *
 121 *     Test the input parameters
 122 *
 123       INFO = 0
 124       UPPER = LSAME( UPLO, 'U' )
 125       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 126          INFO = -1
 127       ELSE IF( N.LT.0 ) THEN
 128          INFO = -2
 129       END IF
 130       IF( INFO.NE.0 ) THEN
 131          CALL XERBLA( 'ZHPTRD', -INFO )
 132          RETURN
 133       END IF
 134 *
 135 *     Quick return if possible
 136 *
 137       IF( N.LE.0 )
 138      $   RETURN
 139 *
 140       IF( UPPER ) THEN
 141 *
 142 *        Reduce the upper triangle of A.
 143 *        I1 is the index in AP of A(1,I+1).
 144 *
 145          I1 = N*( N-1 ) / 2 + 1
 146          AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) )
 147          DO 10 I = N - 1, 1, -1
 148 *
 149 *           Generate elementary reflector H(i) = I - tau * v * v'
 150 *           to annihilate A(1:i-1,i+1)
 151 *
 152             ALPHA = AP( I1+I-1 )
 153             CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI )
 154             E( I ) = ALPHA
 155 *
 156             IF( TAUI.NE.ZERO ) THEN
 157 *
 158 *              Apply H(i) from both sides to A(1:i,1:i)
 159 *
 160                AP( I1+I-1 ) = ONE
 161 *
 162 *              Compute  y := tau * A * v  storing y in TAU(1:i)
 163 *
 164                CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
 165      $                     1 )
 166 *
 167 *              Compute  w := y - 1/2 * tau * (y'*v) * v
 168 *
 169                ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 )
 170                CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
 171 *
 172 *              Apply the transformation as a rank-2 update:
 173 *                 A := A - v * w' - w * v'
 174 *
 175                CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
 176 *
 177             END IF
 178             AP( I1+I-1 ) = E( I )
 179             D( I+1 ) = AP( I1+I )
 180             TAU( I ) = TAUI
 181             I1 = I1 - I
 182    10    CONTINUE
 183          D( 1 ) = AP( 1 )
 184       ELSE
 185 *
 186 *        Reduce the lower triangle of A. II is the index in AP of
 187 *        A(i,i) and I1I1 is the index of A(i+1,i+1).
 188 *
 189          II = 1
 190          AP( 1 ) = DBLE( AP( 1 ) )
 191          DO 20 I = 1, N - 1
 192             I1I1 = II + N - I + 1
 193 *
 194 *           Generate elementary reflector H(i) = I - tau * v * v'
 195 *           to annihilate A(i+2:n,i)
 196 *
 197             ALPHA = AP( II+1 )
 198             CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI )
 199             E( I ) = ALPHA
 200 *
 201             IF( TAUI.NE.ZERO ) THEN
 202 *
 203 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
 204 *
 205                AP( II+1 ) = ONE
 206 *
 207 *              Compute  y := tau * A * v  storing y in TAU(i:n-1)
 208 *
 209                CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
 210      $                     ZERO, TAU( I ), 1 )
 211 *
 212 *              Compute  w := y - 1/2 * tau * (y'*v) * v
 213 *
 214                ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ),
 215      $                 1 )
 216                CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
 217 *
 218 *              Apply the transformation as a rank-2 update:
 219 *                 A := A - v * w' - w * v'
 220 *
 221                CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
 222      $                     AP( I1I1 ) )
 223 *
 224             END IF
 225             AP( II+1 ) = E( I )
 226             D( I ) = AP( II )
 227             TAU( I ) = TAUI
 228             II = I1I1
 229    20    CONTINUE
 230          D( N ) = AP( II )
 231       END IF
 232 *
 233       RETURN
 234 *
 235 *     End of ZHPTRD
 236 *
 237       END