src/LAPACK/dsytrd.f

   1       SUBROUTINE DSYTRD( 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   A( LDA, * ), D( * ), E( * ), TAU( * ),
  13      $                   WORK( * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  DSYTRD reduces a real symmetric matrix A to real symmetric
  20 *  tridiagonal form T by an orthogonal similarity transformation:
  21 *  Q**T * 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) DOUBLE PRECISION array, dimension (LDA,N)
  34 *          On entry, the symmetric 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 orthogonal
  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 orthogonal 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) DOUBLE PRECISION array, dimension (N-1)
  64 *          The scalar factors of the elementary reflectors (see Further
  65 *          Details).
  66 *
  67 *  WORK    (workspace/output) DOUBLE PRECISION 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 real scalar, and v is a real 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 real scalar, and v is a real 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 *     ..
 133 *     .. Local Scalars ..
 134       LOGICAL            LQUERY, UPPER
 135       INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
 136      $                   NBMIN, NX
 137 *     ..
 138 *     .. External Subroutines ..
 139       EXTERNAL           DLATRD, DSYR2K, DSYTD2, XERBLA
 140 *     ..
 141 *     .. Intrinsic Functions ..
 142       INTRINSIC          MAX
 143 *     ..
 144 *     .. External Functions ..
 145       LOGICAL            LSAME
 146       INTEGER            ILAENV
 147       EXTERNAL           LSAME, ILAENV
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Test the input parameters
 152 *
 153       INFO = 0
 154       UPPER = LSAME( UPLO, 'U' )
 155       LQUERY = ( LWORK.EQ.-1 )
 156       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 157          INFO = -1
 158       ELSE IF( N.LT.0 ) THEN
 159          INFO = -2
 160       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 161          INFO = -4
 162       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
 163          INFO = -9
 164       END IF
 165 *
 166       IF( INFO.EQ.0 ) THEN
 167 *
 168 *        Determine the block size.
 169 *
 170          NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
 171          LWKOPT = N*NB
 172          WORK( 1 ) = LWKOPT
 173       END IF
 174 *
 175       IF( INFO.NE.0 ) THEN
 176          CALL XERBLA( 'DSYTRD', -INFO )
 177          RETURN
 178       ELSE IF( LQUERY ) THEN
 179          RETURN
 180       END IF
 181 *
 182 *     Quick return if possible
 183 *
 184       IF( N.EQ.0 ) THEN
 185          WORK( 1 ) = 1
 186          RETURN
 187       END IF
 188 *
 189       NX = N
 190       IWS = 1
 191       IF( NB.GT.1 .AND. NB.LT.N ) THEN
 192 *
 193 *        Determine when to cross over from blocked to unblocked code
 194 *        (last block is always handled by unblocked code).
 195 *
 196          NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
 197          IF( NX.LT.N ) THEN
 198 *
 199 *           Determine if workspace is large enough for blocked code.
 200 *
 201             LDWORK = N
 202             IWS = LDWORK*NB
 203             IF( LWORK.LT.IWS ) THEN
 204 *
 205 *              Not enough workspace to use optimal NB:  determine the
 206 *              minimum value of NB, and reduce NB or force use of
 207 *              unblocked code by setting NX = N.
 208 *
 209                NB = MAX( LWORK / LDWORK, 1 )
 210                NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 )
 211                IF( NB.LT.NBMIN )
 212      $            NX = N
 213             END IF
 214          ELSE
 215             NX = N
 216          END IF
 217       ELSE
 218          NB = 1
 219       END IF
 220 *
 221       IF( UPPER ) THEN
 222 *
 223 *        Reduce the upper triangle of A.
 224 *        Columns 1:kk are handled by the unblocked method.
 225 *
 226          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
 227          DO 20 I = N - NB + 1, KK + 1, -NB
 228 *
 229 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 230 *           matrix W which is needed to update the unreduced part of
 231 *           the matrix
 232 *
 233             CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
 234      $                   LDWORK )
 235 *
 236 *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
 237 *           update of the form:  A := A - V*W' - W*V'
 238 *
 239             CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ),
 240      $                   LDA, WORK, LDWORK, ONE, A, LDA )
 241 *
 242 *           Copy superdiagonal elements back into A, and diagonal
 243 *           elements into D
 244 *
 245             DO 10 J = I, I + NB - 1
 246                A( J-1, J ) = E( J-1 )
 247                D( J ) = A( J, J )
 248    10       CONTINUE
 249    20    CONTINUE
 250 *
 251 *        Use unblocked code to reduce the last or only block
 252 *
 253          CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
 254       ELSE
 255 *
 256 *        Reduce the lower triangle of A
 257 *
 258          DO 40 I = 1, N - NX, NB
 259 *
 260 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
 261 *           matrix W which is needed to update the unreduced part of
 262 *           the matrix
 263 *
 264             CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
 265      $                   TAU( I ), WORK, LDWORK )
 266 *
 267 *           Update the unreduced submatrix A(i+ib:n,i+ib:n), using
 268 *           an update of the form:  A := A - V*W' - W*V'
 269 *
 270             CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE,
 271      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
 272      $                   A( I+NB, I+NB ), LDA )
 273 *
 274 *           Copy subdiagonal elements back into A, and diagonal
 275 *           elements into D
 276 *
 277             DO 30 J = I, I + NB - 1
 278                A( J+1, J ) = E( J )
 279                D( J ) = A( J, J )
 280    30       CONTINUE
 281    40    CONTINUE
 282 *
 283 *        Use unblocked code to reduce the last or only block
 284 *
 285          CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
 286      $                TAU( I ), IINFO )
 287       END IF
 288 *
 289       WORK( 1 ) = LWKOPT
 290       RETURN
 291 *
 292 *     End of DSYTRD
 293 *
 294       END