SRC/dsytrd.f

   1 *> \brief \b DSYTRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSYTRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYTRD( 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   A( LDA, * ), D( * ), E( * ), TAU( * ),
  29 *      $                   WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DSYTRD reduces a real symmetric matrix A to real symmetric
  39 *> tridiagonal form T by an orthogonal similarity transformation:
  40 *> Q**T * 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 DOUBLE PRECISION array, dimension (LDA,N)
  62 *>          On entry, the symmetric 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 orthogonal
  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 orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleSYcomputational
 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**T
 159 *>
 160 *>  where tau is a real scalar, and v is a real 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**T
 172 *>
 173 *>  where tau is a real scalar, and v is a real 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 DSYTRD( 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   A( LDA, * ), D( * ), E( * ), TAU( * ),
 206      $                   WORK( * )
 207 *     ..
 208 *
 209 *  =====================================================================
 210 *
 211 *     .. Parameters ..
 212       DOUBLE PRECISION   ONE
 213       PARAMETER          ( ONE = 1.0D+0 )
 214 *     ..
 215 *     .. Local Scalars ..
 216       LOGICAL            LQUERY, UPPER
 217       INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
 218      $                   NBMIN, NX
 219 *     ..
 220 *     .. External Subroutines ..
 221       EXTERNAL           DLATRD, DSYR2K, DSYTD2, XERBLA
 222 *     ..
 223 *     .. Intrinsic Functions ..
 224       INTRINSIC          MAX
 225 *     ..
 226 *     .. External Functions ..
 227       LOGICAL            LSAME
 228       INTEGER            ILAENV
 229       EXTERNAL           LSAME, ILAENV
 230 *     ..
 231 *     .. Executable Statements ..
 232 *
 233 *     Test the input parameters
 234 *
 235       INFO = 0
 236       UPPER = LSAME( UPLO, 'U' )
 237       LQUERY = ( LWORK.EQ.-1 )
 238       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 239          INFO = -1
 240       ELSE IF( N.LT.0 ) THEN
 241          INFO = -2
 242       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 243          INFO = -4
 244       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
 245          INFO = -9
 246       END IF
 247 *
 248       IF( INFO.EQ.0 ) THEN
 249 *
 250 *        Determine the block size.
 251 *
 252          NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
 253          LWKOPT = N*NB
 254          WORK( 1 ) = LWKOPT
 255       END IF
 256 *
 257       IF( INFO.NE.0 ) THEN
 258          CALL XERBLA( 'DSYTRD', -INFO )
 259          RETURN
 260       ELSE IF( LQUERY ) THEN
 261          RETURN
 262       END IF
 263 *
 264 *     Quick return if possible
 265 *
 266       IF( N.EQ.0 ) THEN
 267          WORK( 1 ) = 1
 268          RETURN
 269       END IF
 270 *
 271       NX = N
 272       IWS = 1
 273       IF( NB.GT.1 .AND. NB.LT.N ) THEN
 274 *
 275 *        Determine when to cross over from blocked to unblocked code
 276 *        (last block is always handled by unblocked code).
 277 *
 278          NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
 279          IF( NX.LT.N ) THEN
 280 *
 281 *           Determine if workspace is large enough for blocked code.
 282 *
 283             LDWORK = N
 284             IWS = LDWORK*NB
 285             IF( LWORK.LT.IWS ) THEN
 286 *
 287 *              Not enough workspace to use optimal NB:  determine the
 288 *              minimum value of NB, and reduce NB or force use of
 289 *              unblocked code by setting NX = N.
 290 *
 291                NB = MAX( LWORK / LDWORK, 1 )
 292                NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 )
 293                IF( NB.LT.NBMIN )
 294      $            NX = N
 295             END IF
 296          ELSE
 297             NX = N
 298          END IF
 299       ELSE
 300          NB = 1
 301       END IF
 302 *
 303       IF( UPPER ) THEN
 304 *
 305 *        Reduce the upper triangle of A.
 306 *        Columns 1:kk are handled by the unblocked method.
 307 *
 308          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
 309          DO 20 I = N - NB + 1, KK + 1, -NB
 310 *
 311 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 312 *           matrix W which is needed to update the unreduced part of
 313 *           the matrix
 314 *
 315             CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
 316      $                   LDWORK )
 317 *
 318 *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
 319 *           update of the form:  A := A - V*W**T - W*V**T
 320 *
 321             CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ),
 322      $                   LDA, WORK, LDWORK, ONE, A, LDA )
 323 *
 324 *           Copy superdiagonal elements back into A, and diagonal
 325 *           elements into D
 326 *
 327             DO 10 J = I, I + NB - 1
 328                A( J-1, J ) = E( J-1 )
 329                D( J ) = A( J, J )
 330    10       CONTINUE
 331    20    CONTINUE
 332 *
 333 *        Use unblocked code to reduce the last or only block
 334 *
 335          CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
 336       ELSE
 337 *
 338 *        Reduce the lower triangle of A
 339 *
 340          DO 40 I = 1, N - NX, NB
 341 *
 342 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 343 *           matrix W which is needed to update the unreduced part of
 344 *           the matrix
 345 *
 346             CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
 347      $                   TAU( I ), WORK, LDWORK )
 348 *
 349 *           Update the unreduced submatrix A(i+ib:n,i+ib:n), using
 350 *           an update of the form:  A := A - V*W**T - W*V**T
 351 *
 352             CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE,
 353      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
 354      $                   A( I+NB, I+NB ), LDA )
 355 *
 356 *           Copy subdiagonal elements back into A, and diagonal
 357 *           elements into D
 358 *
 359             DO 30 J = I, I + NB - 1
 360                A( J+1, J ) = E( J )
 361                D( J ) = A( J, J )
 362    30       CONTINUE
 363    40    CONTINUE
 364 *
 365 *        Use unblocked code to reduce the last or only block
 366 *
 367          CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
 368      $                TAU( I ), IINFO )
 369       END IF
 370 *
 371       WORK( 1 ) = LWKOPT
 372       RETURN
 373 *
 374 *     End of DSYTRD
 375 *
 376       END