001:       SUBROUTINE SLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          UPLO
009:       INTEGER            LDA, LDW, N, NB
010: *     ..
011: *     .. Array Arguments ..
012:       REAL               A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  SLATRD reduces NB rows and columns of a real symmetric matrix A to
019: *  symmetric tridiagonal form by an orthogonal similarity
020: *  transformation Q' * A * Q, and returns the matrices V and W which are
021: *  needed to apply the transformation to the unreduced part of A.
022: *
023: *  If UPLO = 'U', SLATRD reduces the last NB rows and columns of a
024: *  matrix, of which the upper triangle is supplied;
025: *  if UPLO = 'L', SLATRD reduces the first NB rows and columns of a
026: *  matrix, of which the lower triangle is supplied.
027: *
028: *  This is an auxiliary routine called by SSYTRD.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  UPLO    (input) CHARACTER*1
034: *          Specifies whether the upper or lower triangular part of the
035: *          symmetric matrix A is stored:
036: *          = 'U': Upper triangular
037: *          = 'L': Lower triangular
038: *
039: *  N       (input) INTEGER
040: *          The order of the matrix A.
041: *
042: *  NB      (input) INTEGER
043: *          The number of rows and columns to be reduced.
044: *
045: *  A       (input/output) REAL array, dimension (LDA,N)
046: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
047: *          n-by-n upper triangular part of A contains the upper
048: *          triangular part of the matrix A, and the strictly lower
049: *          triangular part of A is not referenced.  If UPLO = 'L', the
050: *          leading n-by-n lower triangular part of A contains the lower
051: *          triangular part of the matrix A, and the strictly upper
052: *          triangular part of A is not referenced.
053: *          On exit:
054: *          if UPLO = 'U', the last NB columns have been reduced to
055: *            tridiagonal form, with the diagonal elements overwriting
056: *            the diagonal elements of A; the elements above the diagonal
057: *            with the array TAU, represent the orthogonal matrix Q as a
058: *            product of elementary reflectors;
059: *          if UPLO = 'L', the first NB columns have been reduced to
060: *            tridiagonal form, with the diagonal elements overwriting
061: *            the diagonal elements of A; the elements below the diagonal
062: *            with the array TAU, represent the  orthogonal matrix Q as a
063: *            product of elementary reflectors.
064: *          See Further Details.
065: *
066: *  LDA     (input) INTEGER
067: *          The leading dimension of the array A.  LDA >= (1,N).
068: *
069: *  E       (output) REAL array, dimension (N-1)
070: *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
071: *          elements of the last NB columns of the reduced matrix;
072: *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
073: *          the first NB columns of the reduced matrix.
074: *
075: *  TAU     (output) REAL array, dimension (N-1)
076: *          The scalar factors of the elementary reflectors, stored in
077: *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
078: *          See Further Details.
079: *
080: *  W       (output) REAL array, dimension (LDW,NB)
081: *          The n-by-nb matrix W required to update the unreduced part
082: *          of A.
083: *
084: *  LDW     (input) INTEGER
085: *          The leading dimension of the array W. LDW >= max(1,N).
086: *
087: *  Further Details
088: *  ===============
089: *
090: *  If UPLO = 'U', the matrix Q is represented as a product of elementary
091: *  reflectors
092: *
093: *     Q = H(n) H(n-1) . . . H(n-nb+1).
094: *
095: *  Each H(i) has the form
096: *
097: *     H(i) = I - tau * v * v'
098: *
099: *  where tau is a real scalar, and v is a real vector with
100: *  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
101: *  and tau in TAU(i-1).
102: *
103: *  If UPLO = 'L', the matrix Q is represented as a product of elementary
104: *  reflectors
105: *
106: *     Q = H(1) H(2) . . . H(nb).
107: *
108: *  Each H(i) has the form
109: *
110: *     H(i) = I - tau * v * v'
111: *
112: *  where tau is a real scalar, and v is a real vector with
113: *  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
114: *  and tau in TAU(i).
115: *
116: *  The elements of the vectors v together form the n-by-nb matrix V
117: *  which is needed, with W, to apply the transformation to the unreduced
118: *  part of the matrix, using a symmetric rank-2k update of the form:
119: *  A := A - V*W' - W*V'.
120: *
121: *  The contents of A on exit are illustrated by the following examples
122: *  with n = 5 and nb = 2:
123: *
124: *  if UPLO = 'U':                       if UPLO = 'L':
125: *
126: *    (  a   a   a   v4  v5 )              (  d                  )
127: *    (      a   a   v4  v5 )              (  1   d              )
128: *    (          a   1   v5 )              (  v1  1   a          )
129: *    (              d   1  )              (  v1  v2  a   a      )
130: *    (                  d  )              (  v1  v2  a   a   a  )
131: *
132: *  where d denotes a diagonal element of the reduced matrix, a denotes
133: *  an element of the original matrix that is unchanged, and vi denotes
134: *  an element of the vector defining H(i).
135: *
136: *  =====================================================================
137: *
138: *     .. Parameters ..
139:       REAL               ZERO, ONE, HALF
140:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, HALF = 0.5E+0 )
141: *     ..
142: *     .. Local Scalars ..
143:       INTEGER            I, IW
144:       REAL               ALPHA
145: *     ..
146: *     .. External Subroutines ..
147:       EXTERNAL           SAXPY, SGEMV, SLARFG, SSCAL, SSYMV
148: *     ..
149: *     .. External Functions ..
150:       LOGICAL            LSAME
151:       REAL               SDOT
152:       EXTERNAL           LSAME, SDOT
153: *     ..
154: *     .. Intrinsic Functions ..
155:       INTRINSIC          MIN
156: *     ..
157: *     .. Executable Statements ..
158: *
159: *     Quick return if possible
160: *
161:       IF( N.LE.0 )
162:      $   RETURN
163: *
164:       IF( LSAME( UPLO, 'U' ) ) THEN
165: *
166: *        Reduce last NB columns of upper triangle
167: *
168:          DO 10 I = N, N - NB + 1, -1
169:             IW = I - N + NB
170:             IF( I.LT.N ) THEN
171: *
172: *              Update A(1:i,i)
173: *
174:                CALL SGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
175:      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
176:                CALL SGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
177:      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
178:             END IF
179:             IF( I.GT.1 ) THEN
180: *
181: *              Generate elementary reflector H(i) to annihilate
182: *              A(1:i-2,i)
183: *
184:                CALL SLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) )
185:                E( I-1 ) = A( I-1, I )
186:                A( I-1, I ) = ONE
187: *
188: *              Compute W(1:i-1,i)
189: *
190:                CALL SSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
191:      $                     ZERO, W( 1, IW ), 1 )
192:                IF( I.LT.N ) THEN
193:                   CALL SGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ),
194:      $                        LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
195:                   CALL SGEMV( 'No transpose', I-1, N-I, -ONE,
196:      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
197:      $                        W( 1, IW ), 1 )
198:                   CALL SGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ),
199:      $                        LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
200:                   CALL SGEMV( 'No transpose', I-1, N-I, -ONE,
201:      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
202:      $                        W( 1, IW ), 1 )
203:                END IF
204:                CALL SSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
205:                ALPHA = -HALF*TAU( I-1 )*SDOT( I-1, W( 1, IW ), 1,
206:      $                 A( 1, I ), 1 )
207:                CALL SAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
208:             END IF
209: *
210:    10    CONTINUE
211:       ELSE
212: *
213: *        Reduce first NB columns of lower triangle
214: *
215:          DO 20 I = 1, NB
216: *
217: *           Update A(i:n,i)
218: *
219:             CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
220:      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
221:             CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
222:      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
223:             IF( I.LT.N ) THEN
224: *
225: *              Generate elementary reflector H(i) to annihilate
226: *              A(i+2:n,i)
227: *
228:                CALL SLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
229:      $                      TAU( I ) )
230:                E( I ) = A( I+1, I )
231:                A( I+1, I ) = ONE
232: *
233: *              Compute W(i+1:n,i)
234: *
235:                CALL SSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
236:      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
237:                CALL SGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW,
238:      $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
239:                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
240:      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
241:                CALL SGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
242:      $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
243:                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
244:      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
245:                CALL SSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
246:                ALPHA = -HALF*TAU( I )*SDOT( N-I, W( I+1, I ), 1,
247:      $                 A( I+1, I ), 1 )
248:                CALL SAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
249:             END IF
250: *
251:    20    CONTINUE
252:       END IF
253: *
254:       RETURN
255: *
256: *     End of SLATRD
257: *
258:       END
259: