LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slatrd.f
Go to the documentation of this file.
00001 *> \brief \b SLATRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLATRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            LDA, LDW, N, NB
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> SLATRD reduces NB rows and columns of a real symmetric matrix A to
00038 *> symmetric tridiagonal form by an orthogonal similarity
00039 *> transformation Q**T * A * Q, and returns the matrices V and W which are
00040 *> needed to apply the transformation to the unreduced part of A.
00041 *>
00042 *> If UPLO = 'U', SLATRD reduces the last NB rows and columns of a
00043 *> matrix, of which the upper triangle is supplied;
00044 *> if UPLO = 'L', SLATRD reduces the first NB rows and columns of a
00045 *> matrix, of which the lower triangle is supplied.
00046 *>
00047 *> This is an auxiliary routine called by SSYTRD.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] UPLO
00054 *> \verbatim
00055 *>          UPLO is CHARACTER*1
00056 *>          Specifies whether the upper or lower triangular part of the
00057 *>          symmetric matrix A is stored:
00058 *>          = 'U': Upper triangular
00059 *>          = 'L': Lower triangular
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix A.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] NB
00069 *> \verbatim
00070 *>          NB is INTEGER
00071 *>          The number of rows and columns to be reduced.
00072 *> \endverbatim
00073 *>
00074 *> \param[in,out] A
00075 *> \verbatim
00076 *>          A is REAL array, dimension (LDA,N)
00077 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
00078 *>          n-by-n upper triangular part of A contains the upper
00079 *>          triangular part of the matrix A, and the strictly lower
00080 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00081 *>          leading n-by-n lower triangular part of A contains the lower
00082 *>          triangular part of the matrix A, and the strictly upper
00083 *>          triangular part of A is not referenced.
00084 *>          On exit:
00085 *>          if UPLO = 'U', the last NB columns have been reduced to
00086 *>            tridiagonal form, with the diagonal elements overwriting
00087 *>            the diagonal elements of A; the elements above the diagonal
00088 *>            with the array TAU, represent the orthogonal matrix Q as a
00089 *>            product of elementary reflectors;
00090 *>          if UPLO = 'L', the first NB columns have been reduced to
00091 *>            tridiagonal form, with the diagonal elements overwriting
00092 *>            the diagonal elements of A; the elements below the diagonal
00093 *>            with the array TAU, represent the  orthogonal matrix Q as a
00094 *>            product of elementary reflectors.
00095 *>          See Further Details.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] LDA
00099 *> \verbatim
00100 *>          LDA is INTEGER
00101 *>          The leading dimension of the array A.  LDA >= (1,N).
00102 *> \endverbatim
00103 *>
00104 *> \param[out] E
00105 *> \verbatim
00106 *>          E is REAL array, dimension (N-1)
00107 *>          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
00108 *>          elements of the last NB columns of the reduced matrix;
00109 *>          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
00110 *>          the first NB columns of the reduced matrix.
00111 *> \endverbatim
00112 *>
00113 *> \param[out] TAU
00114 *> \verbatim
00115 *>          TAU is REAL array, dimension (N-1)
00116 *>          The scalar factors of the elementary reflectors, stored in
00117 *>          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
00118 *>          See Further Details.
00119 *> \endverbatim
00120 *>
00121 *> \param[out] W
00122 *> \verbatim
00123 *>          W is REAL array, dimension (LDW,NB)
00124 *>          The n-by-nb matrix W required to update the unreduced part
00125 *>          of A.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] LDW
00129 *> \verbatim
00130 *>          LDW is INTEGER
00131 *>          The leading dimension of the array W. LDW >= max(1,N).
00132 *> \endverbatim
00133 *
00134 *  Authors:
00135 *  ========
00136 *
00137 *> \author Univ. of Tennessee 
00138 *> \author Univ. of California Berkeley 
00139 *> \author Univ. of Colorado Denver 
00140 *> \author NAG Ltd. 
00141 *
00142 *> \date November 2011
00143 *
00144 *> \ingroup doubleOTHERauxiliary
00145 *
00146 *> \par Further Details:
00147 *  =====================
00148 *>
00149 *> \verbatim
00150 *>
00151 *>  If UPLO = 'U', the matrix Q is represented as a product of elementary
00152 *>  reflectors
00153 *>
00154 *>     Q = H(n) H(n-1) . . . H(n-nb+1).
00155 *>
00156 *>  Each H(i) has the form
00157 *>
00158 *>     H(i) = I - tau * v * v**T
00159 *>
00160 *>  where tau is a real scalar, and v is a real vector with
00161 *>  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
00162 *>  and tau in TAU(i-1).
00163 *>
00164 *>  If UPLO = 'L', the matrix Q is represented as a product of elementary
00165 *>  reflectors
00166 *>
00167 *>     Q = H(1) H(2) . . . H(nb).
00168 *>
00169 *>  Each H(i) has the form
00170 *>
00171 *>     H(i) = I - tau * v * v**T
00172 *>
00173 *>  where tau is a real scalar, and v is a real vector with
00174 *>  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
00175 *>  and tau in TAU(i).
00176 *>
00177 *>  The elements of the vectors v together form the n-by-nb matrix V
00178 *>  which is needed, with W, to apply the transformation to the unreduced
00179 *>  part of the matrix, using a symmetric rank-2k update of the form:
00180 *>  A := A - V*W**T - W*V**T.
00181 *>
00182 *>  The contents of A on exit are illustrated by the following examples
00183 *>  with n = 5 and nb = 2:
00184 *>
00185 *>  if UPLO = 'U':                       if UPLO = 'L':
00186 *>
00187 *>    (  a   a   a   v4  v5 )              (  d                  )
00188 *>    (      a   a   v4  v5 )              (  1   d              )
00189 *>    (          a   1   v5 )              (  v1  1   a          )
00190 *>    (              d   1  )              (  v1  v2  a   a      )
00191 *>    (                  d  )              (  v1  v2  a   a   a  )
00192 *>
00193 *>  where d denotes a diagonal element of the reduced matrix, a denotes
00194 *>  an element of the original matrix that is unchanged, and vi denotes
00195 *>  an element of the vector defining H(i).
00196 *> \endverbatim
00197 *>
00198 *  =====================================================================
00199       SUBROUTINE SLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
00200 *
00201 *  -- LAPACK auxiliary routine (version 3.4.0) --
00202 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00203 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00204 *     November 2011
00205 *
00206 *     .. Scalar Arguments ..
00207       CHARACTER          UPLO
00208       INTEGER            LDA, LDW, N, NB
00209 *     ..
00210 *     .. Array Arguments ..
00211       REAL               A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
00212 *     ..
00213 *
00214 *  =====================================================================
00215 *
00216 *     .. Parameters ..
00217       REAL               ZERO, ONE, HALF
00218       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, HALF = 0.5E+0 )
00219 *     ..
00220 *     .. Local Scalars ..
00221       INTEGER            I, IW
00222       REAL               ALPHA
00223 *     ..
00224 *     .. External Subroutines ..
00225       EXTERNAL           SAXPY, SGEMV, SLARFG, SSCAL, SSYMV
00226 *     ..
00227 *     .. External Functions ..
00228       LOGICAL            LSAME
00229       REAL               SDOT
00230       EXTERNAL           LSAME, SDOT
00231 *     ..
00232 *     .. Intrinsic Functions ..
00233       INTRINSIC          MIN
00234 *     ..
00235 *     .. Executable Statements ..
00236 *
00237 *     Quick return if possible
00238 *
00239       IF( N.LE.0 )
00240      $   RETURN
00241 *
00242       IF( LSAME( UPLO, 'U' ) ) THEN
00243 *
00244 *        Reduce last NB columns of upper triangle
00245 *
00246          DO 10 I = N, N - NB + 1, -1
00247             IW = I - N + NB
00248             IF( I.LT.N ) THEN
00249 *
00250 *              Update A(1:i,i)
00251 *
00252                CALL SGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
00253      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
00254                CALL SGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
00255      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
00256             END IF
00257             IF( I.GT.1 ) THEN
00258 *
00259 *              Generate elementary reflector H(i) to annihilate
00260 *              A(1:i-2,i)
00261 *
00262                CALL SLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) )
00263                E( I-1 ) = A( I-1, I )
00264                A( I-1, I ) = ONE
00265 *
00266 *              Compute W(1:i-1,i)
00267 *
00268                CALL SSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
00269      $                     ZERO, W( 1, IW ), 1 )
00270                IF( I.LT.N ) THEN
00271                   CALL SGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ),
00272      $                        LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
00273                   CALL SGEMV( 'No transpose', I-1, N-I, -ONE,
00274      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
00275      $                        W( 1, IW ), 1 )
00276                   CALL SGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ),
00277      $                        LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
00278                   CALL SGEMV( 'No transpose', I-1, N-I, -ONE,
00279      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
00280      $                        W( 1, IW ), 1 )
00281                END IF
00282                CALL SSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
00283                ALPHA = -HALF*TAU( I-1 )*SDOT( I-1, W( 1, IW ), 1,
00284      $                 A( 1, I ), 1 )
00285                CALL SAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
00286             END IF
00287 *
00288    10    CONTINUE
00289       ELSE
00290 *
00291 *        Reduce first NB columns of lower triangle
00292 *
00293          DO 20 I = 1, NB
00294 *
00295 *           Update A(i:n,i)
00296 *
00297             CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
00298      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
00299             CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
00300      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
00301             IF( I.LT.N ) THEN
00302 *
00303 *              Generate elementary reflector H(i) to annihilate
00304 *              A(i+2:n,i)
00305 *
00306                CALL SLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
00307      $                      TAU( I ) )
00308                E( I ) = A( I+1, I )
00309                A( I+1, I ) = ONE
00310 *
00311 *              Compute W(i+1:n,i)
00312 *
00313                CALL SSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
00314      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
00315                CALL SGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW,
00316      $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
00317                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
00318      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
00319                CALL SGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
00320      $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
00321                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
00322      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
00323                CALL SSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
00324                ALPHA = -HALF*TAU( I )*SDOT( N-I, W( I+1, I ), 1,
00325      $                 A( I+1, I ), 1 )
00326                CALL SAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
00327             END IF
00328 *
00329    20    CONTINUE
00330       END IF
00331 *
00332       RETURN
00333 *
00334 *     End of SLATRD
00335 *
00336       END
 All Files Functions