LAPACK 3.3.1
Linear Algebra PACKage

dsytrd.f

Go to the documentation of this file.
00001       SUBROUTINE DSYTRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, LDA, LWORK, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAU( * ),
00014      $                   WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DSYTRD reduces a real symmetric matrix A to real symmetric
00021 *  tridiagonal form T by an orthogonal similarity transformation:
00022 *  Q**T * A * Q = T.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  UPLO    (input) CHARACTER*1
00028 *          = 'U':  Upper triangle of A is stored;
00029 *          = 'L':  Lower triangle of A is stored.
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.  N >= 0.
00033 *
00034 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00035 *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
00036 *          N-by-N upper triangular part of A contains the upper
00037 *          triangular part of the matrix A, and the strictly lower
00038 *          triangular part of A is not referenced.  If UPLO = 'L', the
00039 *          leading N-by-N lower triangular part of A contains the lower
00040 *          triangular part of the matrix A, and the strictly upper
00041 *          triangular part of A is not referenced.
00042 *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
00043 *          of A are overwritten by the corresponding elements of the
00044 *          tridiagonal matrix T, and the elements above the first
00045 *          superdiagonal, with the array TAU, represent the orthogonal
00046 *          matrix Q as a product of elementary reflectors; if UPLO
00047 *          = 'L', the diagonal and first subdiagonal of A are over-
00048 *          written by the corresponding elements of the tridiagonal
00049 *          matrix T, and the elements below the first subdiagonal, with
00050 *          the array TAU, represent the orthogonal matrix Q as a product
00051 *          of elementary reflectors. See Further Details.
00052 *
00053 *  LDA     (input) INTEGER
00054 *          The leading dimension of the array A.  LDA >= max(1,N).
00055 *
00056 *  D       (output) DOUBLE PRECISION array, dimension (N)
00057 *          The diagonal elements of the tridiagonal matrix T:
00058 *          D(i) = A(i,i).
00059 *
00060 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
00061 *          The off-diagonal elements of the tridiagonal matrix T:
00062 *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
00063 *
00064 *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
00065 *          The scalar factors of the elementary reflectors (see Further
00066 *          Details).
00067 *
00068 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00069 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00070 *
00071 *  LWORK   (input) INTEGER
00072 *          The dimension of the array WORK.  LWORK >= 1.
00073 *          For optimum performance LWORK >= N*NB, where NB is the
00074 *          optimal blocksize.
00075 *
00076 *          If LWORK = -1, then a workspace query is assumed; the routine
00077 *          only calculates the optimal size of the WORK array, returns
00078 *          this value as the first entry of the WORK array, and no error
00079 *          message related to LWORK is issued by XERBLA.
00080 *
00081 *  INFO    (output) INTEGER
00082 *          = 0:  successful exit
00083 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00084 *
00085 *  Further Details
00086 *  ===============
00087 *
00088 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
00089 *  reflectors
00090 *
00091 *     Q = H(n-1) . . . H(2) H(1).
00092 *
00093 *  Each H(i) has the form
00094 *
00095 *     H(i) = I - tau * v * v**T
00096 *
00097 *  where tau is a real scalar, and v is a real vector with
00098 *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
00099 *  A(1:i-1,i+1), and tau in TAU(i).
00100 *
00101 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
00102 *  reflectors
00103 *
00104 *     Q = H(1) H(2) . . . H(n-1).
00105 *
00106 *  Each H(i) has the form
00107 *
00108 *     H(i) = I - tau * v * v**T
00109 *
00110 *  where tau is a real scalar, and v is a real vector with
00111 *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
00112 *  and tau in TAU(i).
00113 *
00114 *  The contents of A on exit are illustrated by the following examples
00115 *  with n = 5:
00116 *
00117 *  if UPLO = 'U':                       if UPLO = 'L':
00118 *
00119 *    (  d   e   v2  v3  v4 )              (  d                  )
00120 *    (      d   e   v3  v4 )              (  e   d              )
00121 *    (          d   e   v4 )              (  v1  e   d          )
00122 *    (              d   e  )              (  v1  v2  e   d      )
00123 *    (                  d  )              (  v1  v2  v3  e   d  )
00124 *
00125 *  where d and e denote diagonal and off-diagonal elements of T, and vi
00126 *  denotes an element of the vector defining H(i).
00127 *
00128 *  =====================================================================
00129 *
00130 *     .. Parameters ..
00131       DOUBLE PRECISION   ONE
00132       PARAMETER          ( ONE = 1.0D+0 )
00133 *     ..
00134 *     .. Local Scalars ..
00135       LOGICAL            LQUERY, UPPER
00136       INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
00137      $                   NBMIN, NX
00138 *     ..
00139 *     .. External Subroutines ..
00140       EXTERNAL           DLATRD, DSYR2K, DSYTD2, XERBLA
00141 *     ..
00142 *     .. Intrinsic Functions ..
00143       INTRINSIC          MAX
00144 *     ..
00145 *     .. External Functions ..
00146       LOGICAL            LSAME
00147       INTEGER            ILAENV
00148       EXTERNAL           LSAME, ILAENV
00149 *     ..
00150 *     .. Executable Statements ..
00151 *
00152 *     Test the input parameters
00153 *
00154       INFO = 0
00155       UPPER = LSAME( UPLO, 'U' )
00156       LQUERY = ( LWORK.EQ.-1 )
00157       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00158          INFO = -1
00159       ELSE IF( N.LT.0 ) THEN
00160          INFO = -2
00161       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00162          INFO = -4
00163       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
00164          INFO = -9
00165       END IF
00166 *
00167       IF( INFO.EQ.0 ) THEN
00168 *
00169 *        Determine the block size.
00170 *
00171          NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
00172          LWKOPT = N*NB
00173          WORK( 1 ) = LWKOPT
00174       END IF
00175 *
00176       IF( INFO.NE.0 ) THEN
00177          CALL XERBLA( 'DSYTRD', -INFO )
00178          RETURN
00179       ELSE IF( LQUERY ) THEN
00180          RETURN
00181       END IF
00182 *
00183 *     Quick return if possible
00184 *
00185       IF( N.EQ.0 ) THEN
00186          WORK( 1 ) = 1
00187          RETURN
00188       END IF
00189 *
00190       NX = N
00191       IWS = 1
00192       IF( NB.GT.1 .AND. NB.LT.N ) THEN
00193 *
00194 *        Determine when to cross over from blocked to unblocked code
00195 *        (last block is always handled by unblocked code).
00196 *
00197          NX = MAX( NB, ILAENV( 3, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
00198          IF( NX.LT.N ) THEN
00199 *
00200 *           Determine if workspace is large enough for blocked code.
00201 *
00202             LDWORK = N
00203             IWS = LDWORK*NB
00204             IF( LWORK.LT.IWS ) THEN
00205 *
00206 *              Not enough workspace to use optimal NB:  determine the
00207 *              minimum value of NB, and reduce NB or force use of
00208 *              unblocked code by setting NX = N.
00209 *
00210                NB = MAX( LWORK / LDWORK, 1 )
00211                NBMIN = ILAENV( 2, 'DSYTRD', UPLO, N, -1, -1, -1 )
00212                IF( NB.LT.NBMIN )
00213      $            NX = N
00214             END IF
00215          ELSE
00216             NX = N
00217          END IF
00218       ELSE
00219          NB = 1
00220       END IF
00221 *
00222       IF( UPPER ) THEN
00223 *
00224 *        Reduce the upper triangle of A.
00225 *        Columns 1:kk are handled by the unblocked method.
00226 *
00227          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
00228          DO 20 I = N - NB + 1, KK + 1, -NB
00229 *
00230 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
00231 *           matrix W which is needed to update the unreduced part of
00232 *           the matrix
00233 *
00234             CALL DLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
00235      $                   LDWORK )
00236 *
00237 *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
00238 *           update of the form:  A := A - V*W**T - W*V**T
00239 *
00240             CALL DSYR2K( UPLO, 'No transpose', I-1, NB, -ONE, A( 1, I ),
00241      $                   LDA, WORK, LDWORK, ONE, A, LDA )
00242 *
00243 *           Copy superdiagonal elements back into A, and diagonal
00244 *           elements into D
00245 *
00246             DO 10 J = I, I + NB - 1
00247                A( J-1, J ) = E( J-1 )
00248                D( J ) = A( J, J )
00249    10       CONTINUE
00250    20    CONTINUE
00251 *
00252 *        Use unblocked code to reduce the last or only block
00253 *
00254          CALL DSYTD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
00255       ELSE
00256 *
00257 *        Reduce the lower triangle of A
00258 *
00259          DO 40 I = 1, N - NX, NB
00260 *
00261 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
00262 *           matrix W which is needed to update the unreduced part of
00263 *           the matrix
00264 *
00265             CALL DLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
00266      $                   TAU( I ), WORK, LDWORK )
00267 *
00268 *           Update the unreduced submatrix A(i+ib:n,i+ib:n), using
00269 *           an update of the form:  A := A - V*W**T - W*V**T
00270 *
00271             CALL DSYR2K( UPLO, 'No transpose', N-I-NB+1, NB, -ONE,
00272      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
00273      $                   A( I+NB, I+NB ), LDA )
00274 *
00275 *           Copy subdiagonal elements back into A, and diagonal
00276 *           elements into D
00277 *
00278             DO 30 J = I, I + NB - 1
00279                A( J+1, J ) = E( J )
00280                D( J ) = A( J, J )
00281    30       CONTINUE
00282    40    CONTINUE
00283 *
00284 *        Use unblocked code to reduce the last or only block
00285 *
00286          CALL DSYTD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
00287      $                TAU( I ), IINFO )
00288       END IF
00289 *
00290       WORK( 1 ) = LWKOPT
00291       RETURN
00292 *
00293 *     End of DSYTRD
00294 *
00295       END
 All Files Functions