src/LAPACK/zhetrd.f

   1       SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, 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, LDA, LWORK, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   D( * ), E( * )
  13       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  ZHETRD reduces a complex Hermitian matrix A to real symmetric
  20 *  tridiagonal form T by a unitary similarity transformation:
  21 *  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 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  34 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
  35 *          N-by-N upper triangular part of A contains the upper
  36 *          triangular part of the matrix A, and the strictly lower
  37 *          triangular part of A is not referenced.  If UPLO = 'L', the
  38 *          leading N-by-N lower triangular part of A contains the lower
  39 *          triangular part of the matrix A, and the strictly upper
  40 *          triangular part of A is not referenced.
  41 *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  42 *          of A are overwritten by the corresponding elements of the
  43 *          tridiagonal matrix T, and the elements above the first
  44 *          superdiagonal, with the array TAU, represent the unitary
  45 *          matrix Q as a product of elementary reflectors; if UPLO
  46 *          = 'L', the diagonal and first subdiagonal of A are over-
  47 *          written by the corresponding elements of the tridiagonal
  48 *          matrix T, and the elements below the first subdiagonal, with
  49 *          the array TAU, represent the unitary matrix Q as a product
  50 *          of elementary reflectors. See Further Details.
  51 *
  52 *  LDA     (input) INTEGER
  53 *          The leading dimension of the array A.  LDA >= max(1,N).
  54 *
  55 *  D       (output) DOUBLE PRECISION array, dimension (N)
  56 *          The diagonal elements of the tridiagonal matrix T:
  57 *          D(i) = A(i,i).
  58 *
  59 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
  60 *          The off-diagonal elements of the tridiagonal matrix T:
  61 *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  62 *
  63 *  TAU     (output) COMPLEX*16 array, dimension (N-1)
  64 *          The scalar factors of the elementary reflectors (see Further
  65 *          Details).
  66 *
  67 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
  68 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  69 *
  70 *  LWORK   (input) INTEGER
  71 *          The dimension of the array WORK.  LWORK >= 1.
  72 *          For optimum performance LWORK >= N*NB, where NB is the
  73 *          optimal blocksize.
  74 *
  75 *          If LWORK = -1, then a workspace query is assumed; the routine
  76 *          only calculates the optimal size of the WORK array, returns
  77 *          this value as the first entry of the WORK array, and no error
  78 *          message related to LWORK is issued by XERBLA.
  79 *
  80 *  INFO    (output) INTEGER
  81 *          = 0:  successful exit
  82 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  83 *
  84 *  Further Details
  85 *  ===============
  86 *
  87 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
  88 *  reflectors
  89 *
  90 *     Q = H(n-1) . . . H(2) H(1).
  91 *
  92 *  Each H(i) has the form
  93 *
  94 *     H(i) = I - tau * v * v'
  95 *
  96 *  where tau is a complex scalar, and v is a complex vector with
  97 *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
  98 *  A(1:i-1,i+1), and tau in TAU(i).
  99 *
 100 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
 101 *  reflectors
 102 *
 103 *     Q = H(1) H(2) . . . H(n-1).
 104 *
 105 *  Each H(i) has the form
 106 *
 107 *     H(i) = I - tau * v * v'
 108 *
 109 *  where tau is a complex scalar, and v is a complex vector with
 110 *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
 111 *  and tau in TAU(i).
 112 *
 113 *  The contents of A on exit are illustrated by the following examples
 114 *  with n = 5:
 115 *
 116 *  if UPLO = 'U':                       if UPLO = 'L':
 117 *
 118 *    (  d   e   v2  v3  v4 )              (  d                  )
 119 *    (      d   e   v3  v4 )              (  e   d              )
 120 *    (          d   e   v4 )              (  v1  e   d          )
 121 *    (              d   e  )              (  v1  v2  e   d      )
 122 *    (                  d  )              (  v1  v2  v3  e   d  )
 123 *
 124 *  where d and e denote diagonal and off-diagonal elements of T, and vi
 125 *  denotes an element of the vector defining H(i).
 126 *
 127 *  =====================================================================
 128 *
 129 *     .. Parameters ..
 130       DOUBLE PRECISION   ONE
 131       PARAMETER          ( ONE = 1.0D+0 )
 132       COMPLEX*16         CONE
 133       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
 134 *     ..
 135 *     .. Local Scalars ..
 136       LOGICAL            LQUERY, UPPER
 137       INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
 138      $                   NBMIN, NX
 139 *     ..
 140 *     .. External Subroutines ..
 141       EXTERNAL           XERBLA, ZHER2K, ZHETD2, ZLATRD
 142 *     ..
 143 *     .. Intrinsic Functions ..
 144       INTRINSIC          MAX
 145 *     ..
 146 *     .. External Functions ..
 147       LOGICAL            LSAME
 148       INTEGER            ILAENV
 149       EXTERNAL           LSAME, ILAENV
 150 *     ..
 151 *     .. Executable Statements ..
 152 *
 153 *     Test the input parameters
 154 *
 155       INFO = 0
 156       UPPER = LSAME( UPLO, 'U' )
 157       LQUERY = ( LWORK.EQ.-1 )
 158       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 159          INFO = -1
 160       ELSE IF( N.LT.0 ) THEN
 161          INFO = -2
 162       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 163          INFO = -4
 164       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
 165          INFO = -9
 166       END IF
 167 *
 168       IF( INFO.EQ.0 ) THEN
 169 *
 170 *        Determine the block size.
 171 *
 172          NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
 173          LWKOPT = N*NB
 174          WORK( 1 ) = LWKOPT
 175       END IF
 176 *
 177       IF( INFO.NE.0 ) THEN
 178          CALL XERBLA( 'ZHETRD', -INFO )
 179          RETURN
 180       ELSE IF( LQUERY ) THEN
 181          RETURN
 182       END IF
 183 *
 184 *     Quick return if possible
 185 *
 186       IF( N.EQ.0 ) THEN
 187          WORK( 1 ) = 1
 188          RETURN
 189       END IF
 190 *
 191       NX = N
 192       IWS = 1
 193       IF( NB.GT.1 .AND. NB.LT.N ) THEN
 194 *
 195 *        Determine when to cross over from blocked to unblocked code
 196 *        (last block is always handled by unblocked code).
 197 *
 198          NX = MAX( NB, ILAENV( 3, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
 199          IF( NX.LT.N ) THEN
 200 *
 201 *           Determine if workspace is large enough for blocked code.
 202 *
 203             LDWORK = N
 204             IWS = LDWORK*NB
 205             IF( LWORK.LT.IWS ) THEN
 206 *
 207 *              Not enough workspace to use optimal NB:  determine the
 208 *              minimum value of NB, and reduce NB or force use of
 209 *              unblocked code by setting NX = N.
 210 *
 211                NB = MAX( LWORK / LDWORK, 1 )
 212                NBMIN = ILAENV( 2, 'ZHETRD', UPLO, N, -1, -1, -1 )
 213                IF( NB.LT.NBMIN )
 214      $            NX = N
 215             END IF
 216          ELSE
 217             NX = N
 218          END IF
 219       ELSE
 220          NB = 1
 221       END IF
 222 *
 223       IF( UPPER ) THEN
 224 *
 225 *        Reduce the upper triangle of A.
 226 *        Columns 1:kk are handled by the unblocked method.
 227 *
 228          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
 229          DO 20 I = N - NB + 1, KK + 1, -NB
 230 *
 231 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 232 *           matrix W which is needed to update the unreduced part of
 233 *           the matrix
 234 *
 235             CALL ZLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
 236      $                   LDWORK )
 237 *
 238 *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
 239 *           update of the form:  A := A - V*W' - W*V'
 240 *
 241             CALL ZHER2K( UPLO, 'No transpose', I-1, NB, -CONE,
 242      $                   A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA )
 243 *
 244 *           Copy superdiagonal elements back into A, and diagonal
 245 *           elements into D
 246 *
 247             DO 10 J = I, I + NB - 1
 248                A( J-1, J ) = E( J-1 )
 249                D( J ) = A( J, J )
 250    10       CONTINUE
 251    20    CONTINUE
 252 *
 253 *        Use unblocked code to reduce the last or only block
 254 *
 255          CALL ZHETD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
 256       ELSE
 257 *
 258 *        Reduce the lower triangle of A
 259 *
 260          DO 40 I = 1, N - NX, NB
 261 *
 262 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 263 *           matrix W which is needed to update the unreduced part of
 264 *           the matrix
 265 *
 266             CALL ZLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
 267      $                   TAU( I ), WORK, LDWORK )
 268 *
 269 *           Update the unreduced submatrix A(i+nb:n,i+nb:n), using
 270 *           an update of the form:  A := A - V*W' - W*V'
 271 *
 272             CALL ZHER2K( UPLO, 'No transpose', N-I-NB+1, NB, -CONE,
 273      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
 274      $                   A( I+NB, I+NB ), LDA )
 275 *
 276 *           Copy subdiagonal elements back into A, and diagonal
 277 *           elements into D
 278 *
 279             DO 30 J = I, I + NB - 1
 280                A( J+1, J ) = E( J )
 281                D( J ) = A( J, J )
 282    30       CONTINUE
 283    40    CONTINUE
 284 *
 285 *        Use unblocked code to reduce the last or only block
 286 *
 287          CALL ZHETD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
 288      $                TAU( I ), IINFO )
 289       END IF
 290 *
 291       WORK( 1 ) = LWKOPT
 292       RETURN
 293 *
 294 *     End of ZHETRD
 295 *
 296       END