LAPACK 3.3.1 Linear Algebra PACKage

# ssytd2.f

Go to the documentation of this file.
```00001       SUBROUTINE SSYTD2( UPLO, N, A, LDA, D, E, TAU, 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, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               A( LDA, * ), D( * ), E( * ), TAU( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  SSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
00020 *  form T by an orthogonal similarity transformation: Q**T * A * Q = T.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  UPLO    (input) CHARACTER*1
00026 *          Specifies whether the upper or lower triangular part of the
00027 *          symmetric matrix A is stored:
00028 *          = 'U':  Upper triangular
00029 *          = 'L':  Lower triangular
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.  N >= 0.
00033 *
00034 *  A       (input/output) REAL 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) REAL array, dimension (N)
00057 *          The diagonal elements of the tridiagonal matrix T:
00058 *          D(i) = A(i,i).
00059 *
00060 *  E       (output) REAL 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) REAL array, dimension (N-1)
00065 *          The scalar factors of the elementary reflectors (see Further
00066 *          Details).
00067 *
00068 *  INFO    (output) INTEGER
00069 *          = 0:  successful exit
00070 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00071 *
00072 *  Further Details
00073 *  ===============
00074 *
00075 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
00076 *  reflectors
00077 *
00078 *     Q = H(n-1) . . . H(2) H(1).
00079 *
00080 *  Each H(i) has the form
00081 *
00082 *     H(i) = I - tau * v * v**T
00083 *
00084 *  where tau is a real scalar, and v is a real vector with
00085 *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
00086 *  A(1:i-1,i+1), and tau in TAU(i).
00087 *
00088 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
00089 *  reflectors
00090 *
00091 *     Q = H(1) H(2) . . . H(n-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(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
00099 *  and tau in TAU(i).
00100 *
00101 *  The contents of A on exit are illustrated by the following examples
00102 *  with n = 5:
00103 *
00104 *  if UPLO = 'U':                       if UPLO = 'L':
00105 *
00106 *    (  d   e   v2  v3  v4 )              (  d                  )
00107 *    (      d   e   v3  v4 )              (  e   d              )
00108 *    (          d   e   v4 )              (  v1  e   d          )
00109 *    (              d   e  )              (  v1  v2  e   d      )
00110 *    (                  d  )              (  v1  v2  v3  e   d  )
00111 *
00112 *  where d and e denote diagonal and off-diagonal elements of T, and vi
00113 *  denotes an element of the vector defining H(i).
00114 *
00115 *  =====================================================================
00116 *
00117 *     .. Parameters ..
00118       REAL               ONE, ZERO, HALF
00119       PARAMETER          ( ONE = 1.0, ZERO = 0.0, HALF = 1.0 / 2.0 )
00120 *     ..
00121 *     .. Local Scalars ..
00122       LOGICAL            UPPER
00123       INTEGER            I
00124       REAL               ALPHA, TAUI
00125 *     ..
00126 *     .. External Subroutines ..
00127       EXTERNAL           SAXPY, SLARFG, SSYMV, SSYR2, XERBLA
00128 *     ..
00129 *     .. External Functions ..
00130       LOGICAL            LSAME
00131       REAL               SDOT
00132       EXTERNAL           LSAME, SDOT
00133 *     ..
00134 *     .. Intrinsic Functions ..
00135       INTRINSIC          MAX, MIN
00136 *     ..
00137 *     .. Executable Statements ..
00138 *
00139 *     Test the input parameters
00140 *
00141       INFO = 0
00142       UPPER = LSAME( UPLO, 'U' )
00143       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00144          INFO = -1
00145       ELSE IF( N.LT.0 ) THEN
00146          INFO = -2
00147       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00148          INFO = -4
00149       END IF
00150       IF( INFO.NE.0 ) THEN
00151          CALL XERBLA( 'SSYTD2', -INFO )
00152          RETURN
00153       END IF
00154 *
00155 *     Quick return if possible
00156 *
00157       IF( N.LE.0 )
00158      \$   RETURN
00159 *
00160       IF( UPPER ) THEN
00161 *
00162 *        Reduce the upper triangle of A
00163 *
00164          DO 10 I = N - 1, 1, -1
00165 *
00166 *           Generate elementary reflector H(i) = I - tau * v * v**T
00167 *           to annihilate A(1:i-1,i+1)
00168 *
00169             CALL SLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI )
00170             E( I ) = A( I, I+1 )
00171 *
00172             IF( TAUI.NE.ZERO ) THEN
00173 *
00174 *              Apply H(i) from both sides to A(1:i,1:i)
00175 *
00176                A( I, I+1 ) = ONE
00177 *
00178 *              Compute  x := tau * A * v  storing x in TAU(1:i)
00179 *
00180                CALL SSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO,
00181      \$                     TAU, 1 )
00182 *
00183 *              Compute  w := x - 1/2 * tau * (x**T * v) * v
00184 *
00185                ALPHA = -HALF*TAUI*SDOT( I, TAU, 1, A( 1, I+1 ), 1 )
00186                CALL SAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 )
00187 *
00188 *              Apply the transformation as a rank-2 update:
00189 *                 A := A - v * w**T - w * v**T
00190 *
00191                CALL SSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A,
00192      \$                     LDA )
00193 *
00194                A( I, I+1 ) = E( I )
00195             END IF
00196             D( I+1 ) = A( I+1, I+1 )
00197             TAU( I ) = TAUI
00198    10    CONTINUE
00199          D( 1 ) = A( 1, 1 )
00200       ELSE
00201 *
00202 *        Reduce the lower triangle of A
00203 *
00204          DO 20 I = 1, N - 1
00205 *
00206 *           Generate elementary reflector H(i) = I - tau * v * v**T
00207 *           to annihilate A(i+2:n,i)
00208 *
00209             CALL SLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
00210      \$                   TAUI )
00211             E( I ) = A( I+1, I )
00212 *
00213             IF( TAUI.NE.ZERO ) THEN
00214 *
00215 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
00216 *
00217                A( I+1, I ) = ONE
00218 *
00219 *              Compute  x := tau * A * v  storing y in TAU(i:n-1)
00220 *
00221                CALL SSYMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA,
00222      \$                     A( I+1, I ), 1, ZERO, TAU( I ), 1 )
00223 *
00224 *              Compute  w := x - 1/2 * tau * (x**T * v) * v
00225 *
00226                ALPHA = -HALF*TAUI*SDOT( N-I, TAU( I ), 1, A( I+1, I ),
00227      \$                 1 )
00228                CALL SAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 )
00229 *
00230 *              Apply the transformation as a rank-2 update:
00231 *                 A := A - v * w**T - w * v**T
00232 *
00233                CALL SSYR2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1,
00234      \$                     A( I+1, I+1 ), LDA )
00235 *
00236                A( I+1, I ) = E( I )
00237             END IF
00238             D( I ) = A( I, I )
00239             TAU( I ) = TAUI
00240    20    CONTINUE
00241          D( N ) = A( N, N )
00242       END IF
00243 *
00244       RETURN
00245 *
00246 *     End of SSYTD2
00247 *
00248       END
```