LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsptrd.f
Go to the documentation of this file.
00001 *> \brief \b DSPTRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DSPTRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSPTRD( UPLO, N, AP, D, E, TAU, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   AP( * ), D( * ), E( * ), TAU( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DSPTRD reduces a real symmetric matrix A stored in packed form to
00038 *> symmetric tridiagonal form T by an orthogonal similarity
00039 *> transformation: Q**T * A * Q = T.
00040 *> \endverbatim
00041 *
00042 *  Arguments:
00043 *  ==========
00044 *
00045 *> \param[in] UPLO
00046 *> \verbatim
00047 *>          UPLO is CHARACTER*1
00048 *>          = 'U':  Upper triangle of A is stored;
00049 *>          = 'L':  Lower triangle of A is stored.
00050 *> \endverbatim
00051 *>
00052 *> \param[in] N
00053 *> \verbatim
00054 *>          N is INTEGER
00055 *>          The order of the matrix A.  N >= 0.
00056 *> \endverbatim
00057 *>
00058 *> \param[in,out] AP
00059 *> \verbatim
00060 *>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
00061 *>          On entry, the upper or lower triangle of the symmetric matrix
00062 *>          A, packed columnwise in a linear array.  The j-th column of A
00063 *>          is stored in the array AP as follows:
00064 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00065 *>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
00066 *>          On exit, if UPLO = 'U', the diagonal and first superdiagonal
00067 *>          of A are overwritten by the corresponding elements of the
00068 *>          tridiagonal matrix T, and the elements above the first
00069 *>          superdiagonal, with the array TAU, represent the orthogonal
00070 *>          matrix Q as a product of elementary reflectors; if UPLO
00071 *>          = 'L', the diagonal and first subdiagonal of A are over-
00072 *>          written by the corresponding elements of the tridiagonal
00073 *>          matrix T, and the elements below the first subdiagonal, with
00074 *>          the array TAU, represent the orthogonal matrix Q as a product
00075 *>          of elementary reflectors. See Further Details.
00076 *> \endverbatim
00077 *>
00078 *> \param[out] D
00079 *> \verbatim
00080 *>          D is DOUBLE PRECISION array, dimension (N)
00081 *>          The diagonal elements of the tridiagonal matrix T:
00082 *>          D(i) = A(i,i).
00083 *> \endverbatim
00084 *>
00085 *> \param[out] E
00086 *> \verbatim
00087 *>          E is DOUBLE PRECISION array, dimension (N-1)
00088 *>          The off-diagonal elements of the tridiagonal matrix T:
00089 *>          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
00090 *> \endverbatim
00091 *>
00092 *> \param[out] TAU
00093 *> \verbatim
00094 *>          TAU is DOUBLE PRECISION array, dimension (N-1)
00095 *>          The scalar factors of the elementary reflectors (see Further
00096 *>          Details).
00097 *> \endverbatim
00098 *>
00099 *> \param[out] INFO
00100 *> \verbatim
00101 *>          INFO is INTEGER
00102 *>          = 0:  successful exit
00103 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00104 *> \endverbatim
00105 *
00106 *  Authors:
00107 *  ========
00108 *
00109 *> \author Univ. of Tennessee 
00110 *> \author Univ. of California Berkeley 
00111 *> \author Univ. of Colorado Denver 
00112 *> \author NAG Ltd. 
00113 *
00114 *> \date November 2011
00115 *
00116 *> \ingroup doubleOTHERcomputational
00117 *
00118 *> \par Further Details:
00119 *  =====================
00120 *>
00121 *> \verbatim
00122 *>
00123 *>  If UPLO = 'U', the matrix Q is represented as a product of elementary
00124 *>  reflectors
00125 *>
00126 *>     Q = H(n-1) . . . H(2) H(1).
00127 *>
00128 *>  Each H(i) has the form
00129 *>
00130 *>     H(i) = I - tau * v * v**T
00131 *>
00132 *>  where tau is a real scalar, and v is a real vector with
00133 *>  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
00134 *>  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
00135 *>
00136 *>  If UPLO = 'L', the matrix Q is represented as a product of elementary
00137 *>  reflectors
00138 *>
00139 *>     Q = H(1) H(2) . . . H(n-1).
00140 *>
00141 *>  Each H(i) has the form
00142 *>
00143 *>     H(i) = I - tau * v * v**T
00144 *>
00145 *>  where tau is a real scalar, and v is a real vector with
00146 *>  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
00147 *>  overwriting A(i+2:n,i), and tau is stored in TAU(i).
00148 *> \endverbatim
00149 *>
00150 *  =====================================================================
00151       SUBROUTINE DSPTRD( UPLO, N, AP, D, E, TAU, INFO )
00152 *
00153 *  -- LAPACK computational routine (version 3.4.0) --
00154 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00155 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00156 *     November 2011
00157 *
00158 *     .. Scalar Arguments ..
00159       CHARACTER          UPLO
00160       INTEGER            INFO, N
00161 *     ..
00162 *     .. Array Arguments ..
00163       DOUBLE PRECISION   AP( * ), D( * ), E( * ), TAU( * )
00164 *     ..
00165 *
00166 *  =====================================================================
00167 *
00168 *     .. Parameters ..
00169       DOUBLE PRECISION   ONE, ZERO, HALF
00170       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0,
00171      $                   HALF = 1.0D0 / 2.0D0 )
00172 *     ..
00173 *     .. Local Scalars ..
00174       LOGICAL            UPPER
00175       INTEGER            I, I1, I1I1, II
00176       DOUBLE PRECISION   ALPHA, TAUI
00177 *     ..
00178 *     .. External Subroutines ..
00179       EXTERNAL           DAXPY, DLARFG, DSPMV, DSPR2, XERBLA
00180 *     ..
00181 *     .. External Functions ..
00182       LOGICAL            LSAME
00183       DOUBLE PRECISION   DDOT
00184       EXTERNAL           LSAME, DDOT
00185 *     ..
00186 *     .. Executable Statements ..
00187 *
00188 *     Test the input parameters
00189 *
00190       INFO = 0
00191       UPPER = LSAME( UPLO, 'U' )
00192       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00193          INFO = -1
00194       ELSE IF( N.LT.0 ) THEN
00195          INFO = -2
00196       END IF
00197       IF( INFO.NE.0 ) THEN
00198          CALL XERBLA( 'DSPTRD', -INFO )
00199          RETURN
00200       END IF
00201 *
00202 *     Quick return if possible
00203 *
00204       IF( N.LE.0 )
00205      $   RETURN
00206 *
00207       IF( UPPER ) THEN
00208 *
00209 *        Reduce the upper triangle of A.
00210 *        I1 is the index in AP of A(1,I+1).
00211 *
00212          I1 = N*( N-1 ) / 2 + 1
00213          DO 10 I = N - 1, 1, -1
00214 *
00215 *           Generate elementary reflector H(i) = I - tau * v * v**T
00216 *           to annihilate A(1:i-1,i+1)
00217 *
00218             CALL DLARFG( I, AP( I1+I-1 ), AP( I1 ), 1, TAUI )
00219             E( I ) = AP( I1+I-1 )
00220 *
00221             IF( TAUI.NE.ZERO ) THEN
00222 *
00223 *              Apply H(i) from both sides to A(1:i,1:i)
00224 *
00225                AP( I1+I-1 ) = ONE
00226 *
00227 *              Compute  y := tau * A * v  storing y in TAU(1:i)
00228 *
00229                CALL DSPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
00230      $                     1 )
00231 *
00232 *              Compute  w := y - 1/2 * tau * (y**T *v) * v
00233 *
00234                ALPHA = -HALF*TAUI*DDOT( I, TAU, 1, AP( I1 ), 1 )
00235                CALL DAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
00236 *
00237 *              Apply the transformation as a rank-2 update:
00238 *                 A := A - v * w**T - w * v**T
00239 *
00240                CALL DSPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
00241 *
00242                AP( I1+I-1 ) = E( I )
00243             END IF
00244             D( I+1 ) = AP( I1+I )
00245             TAU( I ) = TAUI
00246             I1 = I1 - I
00247    10    CONTINUE
00248          D( 1 ) = AP( 1 )
00249       ELSE
00250 *
00251 *        Reduce the lower triangle of A. II is the index in AP of
00252 *        A(i,i) and I1I1 is the index of A(i+1,i+1).
00253 *
00254          II = 1
00255          DO 20 I = 1, N - 1
00256             I1I1 = II + N - I + 1
00257 *
00258 *           Generate elementary reflector H(i) = I - tau * v * v**T
00259 *           to annihilate A(i+2:n,i)
00260 *
00261             CALL DLARFG( N-I, AP( II+1 ), AP( II+2 ), 1, TAUI )
00262             E( I ) = AP( II+1 )
00263 *
00264             IF( TAUI.NE.ZERO ) THEN
00265 *
00266 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
00267 *
00268                AP( II+1 ) = ONE
00269 *
00270 *              Compute  y := tau * A * v  storing y in TAU(i:n-1)
00271 *
00272                CALL DSPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
00273      $                     ZERO, TAU( I ), 1 )
00274 *
00275 *              Compute  w := y - 1/2 * tau * (y**T *v) * v
00276 *
00277                ALPHA = -HALF*TAUI*DDOT( N-I, TAU( I ), 1, AP( II+1 ),
00278      $                 1 )
00279                CALL DAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
00280 *
00281 *              Apply the transformation as a rank-2 update:
00282 *                 A := A - v * w**T - w * v**T
00283 *
00284                CALL DSPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
00285      $                     AP( I1I1 ) )
00286 *
00287                AP( II+1 ) = E( I )
00288             END IF
00289             D( I ) = AP( II )
00290             TAU( I ) = TAUI
00291             II = I1I1
00292    20    CONTINUE
00293          D( N ) = AP( II )
00294       END IF
00295 *
00296       RETURN
00297 *
00298 *     End of DSPTRD
00299 *
00300       END
 All Files Functions