1*> \brief \b SLATRD reduces the first nb rows and columns of a symmetric/Hermitian matrix A to real tridiagonal form by an orthogonal similarity transformation.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLATRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            LDA, LDW, N, NB
26*       ..
27*       .. Array Arguments ..
28*       REAL               A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> SLATRD reduces NB rows and columns of a real symmetric matrix A to
38*> symmetric tridiagonal form by an orthogonal similarity
39*> transformation Q**T * A * Q, and returns the matrices V and W which are
40*> needed to apply the transformation to the unreduced part of A.
41*>
42*> If UPLO = 'U', SLATRD reduces the last NB rows and columns of a
43*> matrix, of which the upper triangle is supplied;
44*> if UPLO = 'L', SLATRD reduces the first NB rows and columns of a
45*> matrix, of which the lower triangle is supplied.
46*>
47*> This is an auxiliary routine called by SSYTRD.
48*> \endverbatim
49*
50*  Arguments:
51*  ==========
52*
53*> \param[in] UPLO
54*> \verbatim
55*>          UPLO is CHARACTER*1
56*>          Specifies whether the upper or lower triangular part of the
57*>          symmetric matrix A is stored:
58*>          = 'U': Upper triangular
59*>          = 'L': Lower triangular
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*>          N is INTEGER
65*>          The order of the matrix A.
66*> \endverbatim
67*>
68*> \param[in] NB
69*> \verbatim
70*>          NB is INTEGER
71*>          The number of rows and columns to be reduced.
72*> \endverbatim
73*>
74*> \param[in,out] A
75*> \verbatim
76*>          A is REAL array, dimension (LDA,N)
77*>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
78*>          n-by-n upper triangular part of A contains the upper
79*>          triangular part of the matrix A, and the strictly lower
80*>          triangular part of A is not referenced.  If UPLO = 'L', the
81*>          leading n-by-n lower triangular part of A contains the lower
82*>          triangular part of the matrix A, and the strictly upper
83*>          triangular part of A is not referenced.
84*>          On exit:
85*>          if UPLO = 'U', the last NB columns have been reduced to
86*>            tridiagonal form, with the diagonal elements overwriting
87*>            the diagonal elements of A; the elements above the diagonal
88*>            with the array TAU, represent the orthogonal matrix Q as a
89*>            product of elementary reflectors;
90*>          if UPLO = 'L', the first NB columns have been reduced to
91*>            tridiagonal form, with the diagonal elements overwriting
92*>            the diagonal elements of A; the elements below the diagonal
93*>            with the array TAU, represent the  orthogonal matrix Q as a
94*>            product of elementary reflectors.
95*>          See Further Details.
96*> \endverbatim
97*>
98*> \param[in] LDA
99*> \verbatim
100*>          LDA is INTEGER
101*>          The leading dimension of the array A.  LDA >= (1,N).
102*> \endverbatim
103*>
104*> \param[out] E
105*> \verbatim
106*>          E is REAL array, dimension (N-1)
107*>          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
108*>          elements of the last NB columns of the reduced matrix;
109*>          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
110*>          the first NB columns of the reduced matrix.
111*> \endverbatim
112*>
113*> \param[out] TAU
114*> \verbatim
115*>          TAU is REAL array, dimension (N-1)
116*>          The scalar factors of the elementary reflectors, stored in
117*>          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
118*>          See Further Details.
119*> \endverbatim
120*>
121*> \param[out] W
122*> \verbatim
123*>          W is REAL array, dimension (LDW,NB)
124*>          The n-by-nb matrix W required to update the unreduced part
125*>          of A.
126*> \endverbatim
127*>
128*> \param[in] LDW
129*> \verbatim
130*>          LDW is INTEGER
131*>          The leading dimension of the array W. LDW >= max(1,N).
132*> \endverbatim
133*
134*  Authors:
135*  ========
136*
137*> \author Univ. of Tennessee
138*> \author Univ. of California Berkeley
139*> \author Univ. of Colorado Denver
140*> \author NAG Ltd.
141*
142*> \date September 2012
143*
144*> \ingroup doubleOTHERauxiliary
145*
146*> \par Further Details:
147*  =====================
148*>
149*> \verbatim
150*>
151*>  If UPLO = 'U', the matrix Q is represented as a product of elementary
152*>  reflectors
153*>
154*>     Q = H(n) H(n-1) . . . H(n-nb+1).
155*>
156*>  Each H(i) has the form
157*>
158*>     H(i) = I - tau * v * v**T
159*>
160*>  where tau is a real scalar, and v is a real vector with
161*>  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
162*>  and tau in TAU(i-1).
163*>
164*>  If UPLO = 'L', the matrix Q is represented as a product of elementary
165*>  reflectors
166*>
167*>     Q = H(1) H(2) . . . H(nb).
168*>
169*>  Each H(i) has the form
170*>
171*>     H(i) = I - tau * v * v**T
172*>
173*>  where tau is a real scalar, and v is a real vector with
174*>  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
175*>  and tau in TAU(i).
176*>
177*>  The elements of the vectors v together form the n-by-nb matrix V
178*>  which is needed, with W, to apply the transformation to the unreduced
179*>  part of the matrix, using a symmetric rank-2k update of the form:
180*>  A := A - V*W**T - W*V**T.
181*>
182*>  The contents of A on exit are illustrated by the following examples
183*>  with n = 5 and nb = 2:
184*>
185*>  if UPLO = 'U':                       if UPLO = 'L':
186*>
187*>    (  a   a   a   v4  v5 )              (  d                  )
188*>    (      a   a   v4  v5 )              (  1   d              )
189*>    (          a   1   v5 )              (  v1  1   a          )
190*>    (              d   1  )              (  v1  v2  a   a      )
191*>    (                  d  )              (  v1  v2  a   a   a  )
192*>
193*>  where d denotes a diagonal element of the reduced matrix, a denotes
194*>  an element of the original matrix that is unchanged, and vi denotes
195*>  an element of the vector defining H(i).
196*> \endverbatim
197*>
198*  =====================================================================
199      SUBROUTINE SLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
200*
201*  -- LAPACK auxiliary routine (version 3.4.2) --
202*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
203*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
204*     September 2012
205*
206*     .. Scalar Arguments ..
207      CHARACTER          UPLO
208      INTEGER            LDA, LDW, N, NB
209*     ..
210*     .. Array Arguments ..
211      REAL               A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
212*     ..
213*
214*  =====================================================================
215*
216*     .. Parameters ..
217      REAL               ZERO, ONE, HALF
218      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, HALF = 0.5E+0 )
219*     ..
220*     .. Local Scalars ..
221      INTEGER            I, IW
222      REAL               ALPHA
223*     ..
224*     .. External Subroutines ..
225      EXTERNAL           SAXPY, SGEMV, SLARFG, SSCAL, SSYMV
226*     ..
227*     .. External Functions ..
228      LOGICAL            LSAME
229      REAL               SDOT
230      EXTERNAL           LSAME, SDOT
231*     ..
232*     .. Intrinsic Functions ..
233      INTRINSIC          MIN
234*     ..
235*     .. Executable Statements ..
236*
237*     Quick return if possible
238*
239      IF( N.LE.0 )
240     $   RETURN
241*
242      IF( LSAME( UPLO, 'U' ) ) THEN
243*
244*        Reduce last NB columns of upper triangle
245*
246         DO 10 I = N, N - NB + 1, -1
247            IW = I - N + NB
248            IF( I.LT.N ) THEN
249*
250*              Update A(1:i,i)
251*
252               CALL SGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
253     $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
254               CALL SGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
255     $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
256            END IF
257            IF( I.GT.1 ) THEN
258*
259*              Generate elementary reflector H(i) to annihilate
260*              A(1:i-2,i)
261*
262               CALL SLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) )
263               E( I-1 ) = A( I-1, I )
264               A( I-1, I ) = ONE
265*
266*              Compute W(1:i-1,i)
267*
268               CALL SSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
269     $                     ZERO, W( 1, IW ), 1 )
270               IF( I.LT.N ) THEN
271                  CALL SGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ),
272     $                        LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
273                  CALL SGEMV( 'No transpose', I-1, N-I, -ONE,
274     $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
275     $                        W( 1, IW ), 1 )
276                  CALL SGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ),
277     $                        LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
278                  CALL SGEMV( 'No transpose', I-1, N-I, -ONE,
279     $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
280     $                        W( 1, IW ), 1 )
281               END IF
282               CALL SSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
283               ALPHA = -HALF*TAU( I-1 )*SDOT( I-1, W( 1, IW ), 1,
284     $                 A( 1, I ), 1 )
285               CALL SAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
286            END IF
287*
288   10    CONTINUE
289      ELSE
290*
291*        Reduce first NB columns of lower triangle
292*
293         DO 20 I = 1, NB
294*
295*           Update A(i:n,i)
296*
297            CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
298     $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
299            CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
300     $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
301            IF( I.LT.N ) THEN
302*
303*              Generate elementary reflector H(i) to annihilate
304*              A(i+2:n,i)
305*
306               CALL SLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
307     $                      TAU( I ) )
308               E( I ) = A( I+1, I )
309               A( I+1, I ) = ONE
310*
311*              Compute W(i+1:n,i)
312*
313               CALL SSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
314     $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
315               CALL SGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW,
316     $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
317               CALL SGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
318     $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
319               CALL SGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
320     $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
321               CALL SGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
322     $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
323               CALL SSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
324               ALPHA = -HALF*TAU( I )*SDOT( N-I, W( I+1, I ), 1,
325     $                 A( I+1, I ), 1 )
326               CALL SAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
327            END IF
328*
329   20    CONTINUE
330      END IF
331*
332      RETURN
333*
334*     End of SLATRD
335*
336      END
337