SRC/ssptrd.f

   1 *> \brief \b SSPTRD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSPTRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssptrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssptrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssptrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSPTRD( UPLO, N, AP, D, E, TAU, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               AP( * ), D( * ), E( * ), TAU( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SSPTRD 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 REAL 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 REAL 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 REAL 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 REAL 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 realOTHERcomputational
 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 SSPTRD( 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       REAL               AP( * ), D( * ), E( * ), TAU( * )
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       REAL               ONE, ZERO, HALF
 170       PARAMETER          ( ONE = 1.0, ZERO = 0.0, HALF = 1.0 / 2.0 )
 171 *     ..
 172 *     .. Local Scalars ..
 173       LOGICAL            UPPER
 174       INTEGER            I, I1, I1I1, II
 175       REAL               ALPHA, TAUI
 176 *     ..
 177 *     .. External Subroutines ..
 178       EXTERNAL           SAXPY, SLARFG, SSPMV, SSPR2, XERBLA
 179 *     ..
 180 *     .. External Functions ..
 181       LOGICAL            LSAME
 182       REAL               SDOT
 183       EXTERNAL           LSAME, SDOT
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187 *     Test the input parameters
 188 *
 189       INFO = 0
 190       UPPER = LSAME( UPLO, 'U' )
 191       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 192          INFO = -1
 193       ELSE IF( N.LT.0 ) THEN
 194          INFO = -2
 195       END IF
 196       IF( INFO.NE.0 ) THEN
 197          CALL XERBLA( 'SSPTRD', -INFO )
 198          RETURN
 199       END IF
 200 *
 201 *     Quick return if possible
 202 *
 203       IF( N.LE.0 )
 204      $   RETURN
 205 *
 206       IF( UPPER ) THEN
 207 *
 208 *        Reduce the upper triangle of A.
 209 *        I1 is the index in AP of A(1,I+1).
 210 *
 211          I1 = N*( N-1 ) / 2 + 1
 212          DO 10 I = N - 1, 1, -1
 213 *
 214 *           Generate elementary reflector H(i) = I - tau * v * v**T
 215 *           to annihilate A(1:i-1,i+1)
 216 *
 217             CALL SLARFG( I, AP( I1+I-1 ), AP( I1 ), 1, TAUI )
 218             E( I ) = AP( I1+I-1 )
 219 *
 220             IF( TAUI.NE.ZERO ) THEN
 221 *
 222 *              Apply H(i) from both sides to A(1:i,1:i)
 223 *
 224                AP( I1+I-1 ) = ONE
 225 *
 226 *              Compute  y := tau * A * v  storing y in TAU(1:i)
 227 *
 228                CALL SSPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
 229      $                     1 )
 230 *
 231 *              Compute  w := y - 1/2 * tau * (y**T *v) * v
 232 *
 233                ALPHA = -HALF*TAUI*SDOT( I, TAU, 1, AP( I1 ), 1 )
 234                CALL SAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
 235 *
 236 *              Apply the transformation as a rank-2 update:
 237 *                 A := A - v * w**T - w * v**T
 238 *
 239                CALL SSPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
 240 *
 241                AP( I1+I-1 ) = E( I )
 242             END IF
 243             D( I+1 ) = AP( I1+I )
 244             TAU( I ) = TAUI
 245             I1 = I1 - I
 246    10    CONTINUE
 247          D( 1 ) = AP( 1 )
 248       ELSE
 249 *
 250 *        Reduce the lower triangle of A. II is the index in AP of
 251 *        A(i,i) and I1I1 is the index of A(i+1,i+1).
 252 *
 253          II = 1
 254          DO 20 I = 1, N - 1
 255             I1I1 = II + N - I + 1
 256 *
 257 *           Generate elementary reflector H(i) = I - tau * v * v**T
 258 *           to annihilate A(i+2:n,i)
 259 *
 260             CALL SLARFG( N-I, AP( II+1 ), AP( II+2 ), 1, TAUI )
 261             E( I ) = AP( II+1 )
 262 *
 263             IF( TAUI.NE.ZERO ) THEN
 264 *
 265 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
 266 *
 267                AP( II+1 ) = ONE
 268 *
 269 *              Compute  y := tau * A * v  storing y in TAU(i:n-1)
 270 *
 271                CALL SSPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
 272      $                     ZERO, TAU( I ), 1 )
 273 *
 274 *              Compute  w := y - 1/2 * tau * (y**T *v) * v
 275 *
 276                ALPHA = -HALF*TAUI*SDOT( N-I, TAU( I ), 1, AP( II+1 ),
 277      $                 1 )
 278                CALL SAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
 279 *
 280 *              Apply the transformation as a rank-2 update:
 281 *                 A := A - v * w**T - w * v**T
 282 *
 283                CALL SSPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
 284      $                     AP( I1I1 ) )
 285 *
 286                AP( II+1 ) = E( I )
 287             END IF
 288             D( I ) = AP( II )
 289             TAU( I ) = TAUI
 290             II = I1I1
 291    20    CONTINUE
 292          D( N ) = AP( II )
 293       END IF
 294 *
 295       RETURN
 296 *
 297 *     End of SSPTRD
 298 *
 299       END