SRC/zhetrd.f

   1 *> \brief \b ZHETRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHETRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, LDA, LWORK, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   D( * ), E( * )
  29 *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> ZHETRD reduces a complex Hermitian matrix A to real symmetric
  39 *> tridiagonal form T by a unitary similarity transformation:
  40 *> 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] A
  60 *> \verbatim
  61 *>          A is COMPLEX*16 array, dimension (LDA,N)
  62 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
  63 *>          N-by-N upper triangular part of A contains the upper
  64 *>          triangular part of the matrix A, and the strictly lower
  65 *>          triangular part of A is not referenced.  If UPLO = 'L', the
  66 *>          leading N-by-N lower triangular part of A contains the lower
  67 *>          triangular part of the matrix A, and the strictly upper
  68 *>          triangular part of A is not referenced.
  69 *>          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  70 *>          of A are overwritten by the corresponding elements of the
  71 *>          tridiagonal matrix T, and the elements above the first
  72 *>          superdiagonal, with the array TAU, represent the unitary
  73 *>          matrix Q as a product of elementary reflectors; if UPLO
  74 *>          = 'L', the diagonal and first subdiagonal of A are over-
  75 *>          written by the corresponding elements of the tridiagonal
  76 *>          matrix T, and the elements below the first subdiagonal, with
  77 *>          the array TAU, represent the unitary matrix Q as a product
  78 *>          of elementary reflectors. See Further Details.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] LDA
  82 *> \verbatim
  83 *>          LDA is INTEGER
  84 *>          The leading dimension of the array A.  LDA >= max(1,N).
  85 *> \endverbatim
  86 *>
  87 *> \param[out] D
  88 *> \verbatim
  89 *>          D is DOUBLE PRECISION array, dimension (N)
  90 *>          The diagonal elements of the tridiagonal matrix T:
  91 *>          D(i) = A(i,i).
  92 *> \endverbatim
  93 *>
  94 *> \param[out] E
  95 *> \verbatim
  96 *>          E is DOUBLE PRECISION array, dimension (N-1)
  97 *>          The off-diagonal elements of the tridiagonal matrix T:
  98 *>          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] TAU
 102 *> \verbatim
 103 *>          TAU is COMPLEX*16 array, dimension (N-1)
 104 *>          The scalar factors of the elementary reflectors (see Further
 105 *>          Details).
 106 *> \endverbatim
 107 *>
 108 *> \param[out] WORK
 109 *> \verbatim
 110 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 111 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LWORK
 115 *> \verbatim
 116 *>          LWORK is INTEGER
 117 *>          The dimension of the array WORK.  LWORK >= 1.
 118 *>          For optimum performance LWORK >= N*NB, where NB is the
 119 *>          optimal blocksize.
 120 *>
 121 *>          If LWORK = -1, then a workspace query is assumed; the routine
 122 *>          only calculates the optimal size of the WORK array, returns
 123 *>          this value as the first entry of the WORK array, and no error
 124 *>          message related to LWORK is issued by XERBLA.
 125 *> \endverbatim
 126 *>
 127 *> \param[out] INFO
 128 *> \verbatim
 129 *>          INFO is INTEGER
 130 *>          = 0:  successful exit
 131 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 132 *> \endverbatim
 133 *
 134 *  Authors:
 135 *  ========
 136 *
 137 *> \author Univ. of Tennessee
 138 *> \author Univ. of California Berkeley
 139 *> \author Univ. of Colorado Denver
 140 *> \author NAG Ltd.
 141 *
 142 *> \date November 2011
 143 *
 144 *> \ingroup complex16HEcomputational
 145 *
 146 *> \par Further Details:
 147 *  =====================
 148 *>
 149 *> \verbatim
 150 *>
 151 *>  If UPLO = 'U', the matrix Q is represented as a product of elementary
 152 *>  reflectors
 153 *>
 154 *>     Q = H(n-1) . . . H(2) H(1).
 155 *>
 156 *>  Each H(i) has the form
 157 *>
 158 *>     H(i) = I - tau * v * v**H
 159 *>
 160 *>  where tau is a complex scalar, and v is a complex vector with
 161 *>  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
 162 *>  A(1:i-1,i+1), and tau in TAU(i).
 163 *>
 164 *>  If UPLO = 'L', the matrix Q is represented as a product of elementary
 165 *>  reflectors
 166 *>
 167 *>     Q = H(1) H(2) . . . H(n-1).
 168 *>
 169 *>  Each H(i) has the form
 170 *>
 171 *>     H(i) = I - tau * v * v**H
 172 *>
 173 *>  where tau is a complex scalar, and v is a complex vector with
 174 *>  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
 175 *>  and tau in TAU(i).
 176 *>
 177 *>  The contents of A on exit are illustrated by the following examples
 178 *>  with n = 5:
 179 *>
 180 *>  if UPLO = 'U':                       if UPLO = 'L':
 181 *>
 182 *>    (  d   e   v2  v3  v4 )              (  d                  )
 183 *>    (      d   e   v3  v4 )              (  e   d              )
 184 *>    (          d   e   v4 )              (  v1  e   d          )
 185 *>    (              d   e  )              (  v1  v2  e   d      )
 186 *>    (                  d  )              (  v1  v2  v3  e   d  )
 187 *>
 188 *>  where d and e denote diagonal and off-diagonal elements of T, and vi
 189 *>  denotes an element of the vector defining H(i).
 190 *> \endverbatim
 191 *>
 192 *  =====================================================================
 193       SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
 194 *
 195 *  -- LAPACK computational routine (version 3.4.0) --
 196 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 197 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 198 *     November 2011
 199 *
 200 *     .. Scalar Arguments ..
 201       CHARACTER          UPLO
 202       INTEGER            INFO, LDA, LWORK, N
 203 *     ..
 204 *     .. Array Arguments ..
 205       DOUBLE PRECISION   D( * ), E( * )
 206       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
 207 *     ..
 208 *
 209 *  =====================================================================
 210 *
 211 *     .. Parameters ..
 212       DOUBLE PRECISION   ONE
 213       PARAMETER          ( ONE = 1.0D+0 )
 214       COMPLEX*16         CONE
 215       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
 216 *     ..
 217 *     .. Local Scalars ..
 218       LOGICAL            LQUERY, UPPER
 219       INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
 220      $                   NBMIN, NX
 221 *     ..
 222 *     .. External Subroutines ..
 223       EXTERNAL           XERBLA, ZHER2K, ZHETD2, ZLATRD
 224 *     ..
 225 *     .. Intrinsic Functions ..
 226       INTRINSIC          MAX
 227 *     ..
 228 *     .. External Functions ..
 229       LOGICAL            LSAME
 230       INTEGER            ILAENV
 231       EXTERNAL           LSAME, ILAENV
 232 *     ..
 233 *     .. Executable Statements ..
 234 *
 235 *     Test the input parameters
 236 *
 237       INFO = 0
 238       UPPER = LSAME( UPLO, 'U' )
 239       LQUERY = ( LWORK.EQ.-1 )
 240       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 241          INFO = -1
 242       ELSE IF( N.LT.0 ) THEN
 243          INFO = -2
 244       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 245          INFO = -4
 246       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
 247          INFO = -9
 248       END IF
 249 *
 250       IF( INFO.EQ.0 ) THEN
 251 *
 252 *        Determine the block size.
 253 *
 254          NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
 255          LWKOPT = N*NB
 256          WORK( 1 ) = LWKOPT
 257       END IF
 258 *
 259       IF( INFO.NE.0 ) THEN
 260          CALL XERBLA( 'ZHETRD', -INFO )
 261          RETURN
 262       ELSE IF( LQUERY ) THEN
 263          RETURN
 264       END IF
 265 *
 266 *     Quick return if possible
 267 *
 268       IF( N.EQ.0 ) THEN
 269          WORK( 1 ) = 1
 270          RETURN
 271       END IF
 272 *
 273       NX = N
 274       IWS = 1
 275       IF( NB.GT.1 .AND. NB.LT.N ) THEN
 276 *
 277 *        Determine when to cross over from blocked to unblocked code
 278 *        (last block is always handled by unblocked code).
 279 *
 280          NX = MAX( NB, ILAENV( 3, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
 281          IF( NX.LT.N ) THEN
 282 *
 283 *           Determine if workspace is large enough for blocked code.
 284 *
 285             LDWORK = N
 286             IWS = LDWORK*NB
 287             IF( LWORK.LT.IWS ) THEN
 288 *
 289 *              Not enough workspace to use optimal NB:  determine the
 290 *              minimum value of NB, and reduce NB or force use of
 291 *              unblocked code by setting NX = N.
 292 *
 293                NB = MAX( LWORK / LDWORK, 1 )
 294                NBMIN = ILAENV( 2, 'ZHETRD', UPLO, N, -1, -1, -1 )
 295                IF( NB.LT.NBMIN )
 296      $            NX = N
 297             END IF
 298          ELSE
 299             NX = N
 300          END IF
 301       ELSE
 302          NB = 1
 303       END IF
 304 *
 305       IF( UPPER ) THEN
 306 *
 307 *        Reduce the upper triangle of A.
 308 *        Columns 1:kk are handled by the unblocked method.
 309 *
 310          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
 311          DO 20 I = N - NB + 1, KK + 1, -NB
 312 *
 313 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 314 *           matrix W which is needed to update the unreduced part of
 315 *           the matrix
 316 *
 317             CALL ZLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
 318      $                   LDWORK )
 319 *
 320 *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
 321 *           update of the form:  A := A - V*W**H - W*V**H
 322 *
 323             CALL ZHER2K( UPLO, 'No transpose', I-1, NB, -CONE,
 324      $                   A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA )
 325 *
 326 *           Copy superdiagonal elements back into A, and diagonal
 327 *           elements into D
 328 *
 329             DO 10 J = I, I + NB - 1
 330                A( J-1, J ) = E( J-1 )
 331                D( J ) = A( J, J )
 332    10       CONTINUE
 333    20    CONTINUE
 334 *
 335 *        Use unblocked code to reduce the last or only block
 336 *
 337          CALL ZHETD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
 338       ELSE
 339 *
 340 *        Reduce the lower triangle of A
 341 *
 342          DO 40 I = 1, N - NX, NB
 343 *
 344 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 345 *           matrix W which is needed to update the unreduced part of
 346 *           the matrix
 347 *
 348             CALL ZLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
 349      $                   TAU( I ), WORK, LDWORK )
 350 *
 351 *           Update the unreduced submatrix A(i+nb:n,i+nb:n), using
 352 *           an update of the form:  A := A - V*W**H - W*V**H
 353 *
 354             CALL ZHER2K( UPLO, 'No transpose', N-I-NB+1, NB, -CONE,
 355      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
 356      $                   A( I+NB, I+NB ), LDA )
 357 *
 358 *           Copy subdiagonal elements back into A, and diagonal
 359 *           elements into D
 360 *
 361             DO 30 J = I, I + NB - 1
 362                A( J+1, J ) = E( J )
 363                D( J ) = A( J, J )
 364    30       CONTINUE
 365    40    CONTINUE
 366 *
 367 *        Use unblocked code to reduce the last or only block
 368 *
 369          CALL ZHETD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
 370      $                TAU( I ), IINFO )
 371       END IF
 372 *
 373       WORK( 1 ) = LWKOPT
 374       RETURN
 375 *
 376 *     End of ZHETRD
 377 *
 378       END