SRC/dsptrd.f

   1 *> \brief \b DSPTRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSPTRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSPTRD( UPLO, N, AP, D, E, TAU, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   AP( * ), D( * ), E( * ), TAU( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DSPTRD reduces a real symmetric matrix A stored in packed form to
  38 *> symmetric tridiagonal form T by an orthogonal similarity
  39 *> transformation: Q**T * A * Q = T.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] UPLO
  46 *> \verbatim
  47 *>          UPLO is CHARACTER*1
  48 *>          = 'U':  Upper triangle of A is stored;
  49 *>          = 'L':  Lower triangle of A is stored.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] N
  53 *> \verbatim
  54 *>          N is INTEGER
  55 *>          The order of the matrix A.  N >= 0.
  56 *> \endverbatim
  57 *>
  58 *> \param[in,out] AP
  59 *> \verbatim
  60 *>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
  61 *>          On entry, the upper or lower triangle of the symmetric matrix
  62 *>          A, packed columnwise in a linear array.  The j-th column of A
  63 *>          is stored in the array AP as follows:
  64 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  65 *>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
  66 *>          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  67 *>          of A are overwritten by the corresponding elements of the
  68 *>          tridiagonal matrix T, and the elements above the first
  69 *>          superdiagonal, with the array TAU, represent the orthogonal
  70 *>          matrix Q as a product of elementary reflectors; if UPLO
  71 *>          = 'L', the diagonal and first subdiagonal of A are over-
  72 *>          written by the corresponding elements of the tridiagonal
  73 *>          matrix T, and the elements below the first subdiagonal, with
  74 *>          the array TAU, represent the orthogonal matrix Q as a product
  75 *>          of elementary reflectors. See Further Details.
  76 *> \endverbatim
  77 *>
  78 *> \param[out] D
  79 *> \verbatim
  80 *>          D is DOUBLE PRECISION array, dimension (N)
  81 *>          The diagonal elements of the tridiagonal matrix T:
  82 *>          D(i) = A(i,i).
  83 *> \endverbatim
  84 *>
  85 *> \param[out] E
  86 *> \verbatim
  87 *>          E is DOUBLE PRECISION array, dimension (N-1)
  88 *>          The off-diagonal elements of the tridiagonal matrix T:
  89 *>          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  90 *> \endverbatim
  91 *>
  92 *> \param[out] TAU
  93 *> \verbatim
  94 *>          TAU is DOUBLE PRECISION array, dimension (N-1)
  95 *>          The scalar factors of the elementary reflectors (see Further
  96 *>          Details).
  97 *> \endverbatim
  98 *>
  99 *> \param[out] INFO
 100 *> \verbatim
 101 *>          INFO is INTEGER
 102 *>          = 0:  successful exit
 103 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 104 *> \endverbatim
 105 *
 106 *  Authors:
 107 *  ========
 108 *
 109 *> \author Univ. of Tennessee
 110 *> \author Univ. of California Berkeley
 111 *> \author Univ. of Colorado Denver
 112 *> \author NAG Ltd.
 113 *
 114 *> \date November 2011
 115 *
 116 *> \ingroup doubleOTHERcomputational
 117 *
 118 *> \par Further Details:
 119 *  =====================
 120 *>
 121 *> \verbatim
 122 *>
 123 *>  If UPLO = 'U', the matrix Q is represented as a product of elementary
 124 *>  reflectors
 125 *>
 126 *>     Q = H(n-1) . . . H(2) H(1).
 127 *>
 128 *>  Each H(i) has the form
 129 *>
 130 *>     H(i) = I - tau * v * v**T
 131 *>
 132 *>  where tau is a real scalar, and v is a real vector with
 133 *>  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
 134 *>  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
 135 *>
 136 *>  If UPLO = 'L', the matrix Q is represented as a product of elementary
 137 *>  reflectors
 138 *>
 139 *>     Q = H(1) H(2) . . . H(n-1).
 140 *>
 141 *>  Each H(i) has the form
 142 *>
 143 *>     H(i) = I - tau * v * v**T
 144 *>
 145 *>  where tau is a real scalar, and v is a real vector with
 146 *>  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
 147 *>  overwriting A(i+2:n,i), and tau is stored in TAU(i).
 148 *> \endverbatim
 149 *>
 150 *  =====================================================================
 151       SUBROUTINE DSPTRD( UPLO, N, AP, D, E, TAU, INFO )
 152 *
 153 *  -- LAPACK computational routine (version 3.4.0) --
 154 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 155 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 156 *     November 2011
 157 *
 158 *     .. Scalar Arguments ..
 159       CHARACTER          UPLO
 160       INTEGER            INFO, N
 161 *     ..
 162 *     .. Array Arguments ..
 163       DOUBLE PRECISION   AP( * ), D( * ), E( * ), TAU( * )
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       DOUBLE PRECISION   ONE, ZERO, HALF
 170       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0,
 171      $                   HALF = 1.0D0 / 2.0D0 )
 172 *     ..
 173 *     .. Local Scalars ..
 174       LOGICAL            UPPER
 175       INTEGER            I, I1, I1I1, II
 176       DOUBLE PRECISION   ALPHA, TAUI
 177 *     ..
 178 *     .. External Subroutines ..
 179       EXTERNAL           DAXPY, DLARFG, DSPMV, DSPR2, XERBLA
 180 *     ..
 181 *     .. External Functions ..
 182       LOGICAL            LSAME
 183       DOUBLE PRECISION   DDOT
 184       EXTERNAL           LSAME, DDOT
 185 *     ..
 186 *     .. Executable Statements ..
 187 *
 188 *     Test the input parameters
 189 *
 190       INFO = 0
 191       UPPER = LSAME( UPLO, 'U' )
 192       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 193          INFO = -1
 194       ELSE IF( N.LT.0 ) THEN
 195          INFO = -2
 196       END IF
 197       IF( INFO.NE.0 ) THEN
 198          CALL XERBLA( 'DSPTRD', -INFO )
 199          RETURN
 200       END IF
 201 *
 202 *     Quick return if possible
 203 *
 204       IF( N.LE.0 )
 205      $   RETURN
 206 *
 207       IF( UPPER ) THEN
 208 *
 209 *        Reduce the upper triangle of A.
 210 *        I1 is the index in AP of A(1,I+1).
 211 *
 212          I1 = N*( N-1 ) / 2 + 1
 213          DO 10 I = N - 1, 1, -1
 214 *
 215 *           Generate elementary reflector H(i) = I - tau * v * v**T
 216 *           to annihilate A(1:i-1,i+1)
 217 *
 218             CALL DLARFG( I, AP( I1+I-1 ), AP( I1 ), 1, TAUI )
 219             E( I ) = AP( I1+I-1 )
 220 *
 221             IF( TAUI.NE.ZERO ) THEN
 222 *
 223 *              Apply H(i) from both sides to A(1:i,1:i)
 224 *
 225                AP( I1+I-1 ) = ONE
 226 *
 227 *              Compute  y := tau * A * v  storing y in TAU(1:i)
 228 *
 229                CALL DSPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
 230      $                     1 )
 231 *
 232 *              Compute  w := y - 1/2 * tau * (y**T *v) * v
 233 *
 234                ALPHA = -HALF*TAUI*DDOT( I, TAU, 1, AP( I1 ), 1 )
 235                CALL DAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
 236 *
 237 *              Apply the transformation as a rank-2 update:
 238 *                 A := A - v * w**T - w * v**T
 239 *
 240                CALL DSPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
 241 *
 242                AP( I1+I-1 ) = E( I )
 243             END IF
 244             D( I+1 ) = AP( I1+I )
 245             TAU( I ) = TAUI
 246             I1 = I1 - I
 247    10    CONTINUE
 248          D( 1 ) = AP( 1 )
 249       ELSE
 250 *
 251 *        Reduce the lower triangle of A. II is the index in AP of
 252 *        A(i,i) and I1I1 is the index of A(i+1,i+1).
 253 *
 254          II = 1
 255          DO 20 I = 1, N - 1
 256             I1I1 = II + N - I + 1
 257 *
 258 *           Generate elementary reflector H(i) = I - tau * v * v**T
 259 *           to annihilate A(i+2:n,i)
 260 *
 261             CALL DLARFG( N-I, AP( II+1 ), AP( II+2 ), 1, TAUI )
 262             E( I ) = AP( II+1 )
 263 *
 264             IF( TAUI.NE.ZERO ) THEN
 265 *
 266 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
 267 *
 268                AP( II+1 ) = ONE
 269 *
 270 *              Compute  y := tau * A * v  storing y in TAU(i:n-1)
 271 *
 272                CALL DSPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
 273      $                     ZERO, TAU( I ), 1 )
 274 *
 275 *              Compute  w := y - 1/2 * tau * (y**T *v) * v
 276 *
 277                ALPHA = -HALF*TAUI*DDOT( N-I, TAU( I ), 1, AP( II+1 ),
 278      $                 1 )
 279                CALL DAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
 280 *
 281 *              Apply the transformation as a rank-2 update:
 282 *                 A := A - v * w**T - w * v**T
 283 *
 284                CALL DSPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
 285      $                     AP( I1I1 ) )
 286 *
 287                AP( II+1 ) = E( I )
 288             END IF
 289             D( I ) = AP( II )
 290             TAU( I ) = TAUI
 291             II = I1I1
 292    20    CONTINUE
 293          D( N ) = AP( II )
 294       END IF
 295 *
 296       RETURN
 297 *
 298 *     End of DSPTRD
 299 *
 300       END