1: SUBROUTINE DLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER LDA, LDW, N, NB
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * DLATRD reduces NB rows and columns of a real symmetric matrix A to
20: * symmetric tridiagonal form by an orthogonal similarity
21: * transformation Q' * A * Q, and returns the matrices V and W which are
22: * needed to apply the transformation to the unreduced part of A.
23: *
24: * If UPLO = 'U', DLATRD reduces the last NB rows and columns of a
25: * matrix, of which the upper triangle is supplied;
26: * if UPLO = 'L', DLATRD reduces the first NB rows and columns of a
27: * matrix, of which the lower triangle is supplied.
28: *
29: * This is an auxiliary routine called by DSYTRD.
30: *
31: * Arguments
32: * =========
33: *
34: * UPLO (input) CHARACTER*1
35: * Specifies whether the upper or lower triangular part of the
36: * symmetric matrix A is stored:
37: * = 'U': Upper triangular
38: * = 'L': Lower triangular
39: *
40: * N (input) INTEGER
41: * The order of the matrix A.
42: *
43: * NB (input) INTEGER
44: * The number of rows and columns to be reduced.
45: *
46: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
47: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
48: * n-by-n upper triangular part of A contains the upper
49: * triangular part of the matrix A, and the strictly lower
50: * triangular part of A is not referenced. If UPLO = 'L', the
51: * leading n-by-n lower triangular part of A contains the lower
52: * triangular part of the matrix A, and the strictly upper
53: * triangular part of A is not referenced.
54: * On exit:
55: * if UPLO = 'U', the last NB columns have been reduced to
56: * tridiagonal form, with the diagonal elements overwriting
57: * the diagonal elements of A; the elements above the diagonal
58: * with the array TAU, represent the orthogonal matrix Q as a
59: * product of elementary reflectors;
60: * if UPLO = 'L', the first NB columns have been reduced to
61: * tridiagonal form, with the diagonal elements overwriting
62: * the diagonal elements of A; the elements below the diagonal
63: * with the array TAU, represent the orthogonal matrix Q as a
64: * product of elementary reflectors.
65: * See Further Details.
66: *
67: * LDA (input) INTEGER
68: * The leading dimension of the array A. LDA >= (1,N).
69: *
70: * E (output) DOUBLE PRECISION array, dimension (N-1)
71: * If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
72: * elements of the last NB columns of the reduced matrix;
73: * if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
74: * the first NB columns of the reduced matrix.
75: *
76: * TAU (output) DOUBLE PRECISION array, dimension (N-1)
77: * The scalar factors of the elementary reflectors, stored in
78: * TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
79: * See Further Details.
80: *
81: * W (output) DOUBLE PRECISION array, dimension (LDW,NB)
82: * The n-by-nb matrix W required to update the unreduced part
83: * of A.
84: *
85: * LDW (input) INTEGER
86: * The leading dimension of the array W. LDW >= max(1,N).
87: *
88: * Further Details
89: * ===============
90: *
91: * If UPLO = 'U', the matrix Q is represented as a product of elementary
92: * reflectors
93: *
94: * Q = H(n) H(n-1) . . . H(n-nb+1).
95: *
96: * Each H(i) has the form
97: *
98: * H(i) = I - tau * v * v'
99: *
100: * where tau is a real scalar, and v is a real vector with
101: * v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
102: * and tau in TAU(i-1).
103: *
104: * If UPLO = 'L', the matrix Q is represented as a product of elementary
105: * reflectors
106: *
107: * Q = H(1) H(2) . . . H(nb).
108: *
109: * Each H(i) has the form
110: *
111: * H(i) = I - tau * v * v'
112: *
113: * where tau is a real scalar, and v is a real vector with
114: * v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
115: * and tau in TAU(i).
116: *
117: * The elements of the vectors v together form the n-by-nb matrix V
118: * which is needed, with W, to apply the transformation to the unreduced
119: * part of the matrix, using a symmetric rank-2k update of the form:
120: * A := A - V*W' - W*V'.
121: *
122: * The contents of A on exit are illustrated by the following examples
123: * with n = 5 and nb = 2:
124: *
125: * if UPLO = 'U': if UPLO = 'L':
126: *
127: * ( a a a v4 v5 ) ( d )
128: * ( a a v4 v5 ) ( 1 d )
129: * ( a 1 v5 ) ( v1 1 a )
130: * ( d 1 ) ( v1 v2 a a )
131: * ( d ) ( v1 v2 a a a )
132: *
133: * where d denotes a diagonal element of the reduced matrix, a denotes
134: * an element of the original matrix that is unchanged, and vi denotes
135: * an element of the vector defining H(i).
136: *
137: * =====================================================================
138: *
139: * .. Parameters ..
140: DOUBLE PRECISION ZERO, ONE, HALF
141: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
142: * ..
143: * .. Local Scalars ..
144: INTEGER I, IW
145: DOUBLE PRECISION ALPHA
146: * ..
147: * .. External Subroutines ..
148: EXTERNAL DAXPY, DGEMV, DLARFG, DSCAL, DSYMV
149: * ..
150: * .. External Functions ..
151: LOGICAL LSAME
152: DOUBLE PRECISION DDOT
153: EXTERNAL LSAME, DDOT
154: * ..
155: * .. Intrinsic Functions ..
156: INTRINSIC MIN
157: * ..
158: * .. Executable Statements ..
159: *
160: * Quick return if possible
161: *
162: IF( N.LE.0 )
163: $ RETURN
164: *
165: IF( LSAME( UPLO, 'U' ) ) THEN
166: *
167: * Reduce last NB columns of upper triangle
168: *
169: DO 10 I = N, N - NB + 1, -1
170: IW = I - N + NB
171: IF( I.LT.N ) THEN
172: *
173: * Update A(1:i,i)
174: *
175: CALL DGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
176: $ LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
177: CALL DGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
178: $ LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
179: END IF
180: IF( I.GT.1 ) THEN
181: *
182: * Generate elementary reflector H(i) to annihilate
183: * A(1:i-2,i)
184: *
185: CALL DLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) )
186: E( I-1 ) = A( I-1, I )
187: A( I-1, I ) = ONE
188: *
189: * Compute W(1:i-1,i)
190: *
191: CALL DSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
192: $ ZERO, W( 1, IW ), 1 )
193: IF( I.LT.N ) THEN
194: CALL DGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ),
195: $ LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
196: CALL DGEMV( 'No transpose', I-1, N-I, -ONE,
197: $ A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
198: $ W( 1, IW ), 1 )
199: CALL DGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ),
200: $ LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
201: CALL DGEMV( 'No transpose', I-1, N-I, -ONE,
202: $ W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
203: $ W( 1, IW ), 1 )
204: END IF
205: CALL DSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
206: ALPHA = -HALF*TAU( I-1 )*DDOT( I-1, W( 1, IW ), 1,
207: $ A( 1, I ), 1 )
208: CALL DAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
209: END IF
210: *
211: 10 CONTINUE
212: ELSE
213: *
214: * Reduce first NB columns of lower triangle
215: *
216: DO 20 I = 1, NB
217: *
218: * Update A(i:n,i)
219: *
220: CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
221: $ LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
222: CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
223: $ LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
224: IF( I.LT.N ) THEN
225: *
226: * Generate elementary reflector H(i) to annihilate
227: * A(i+2:n,i)
228: *
229: CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
230: $ TAU( I ) )
231: E( I ) = A( I+1, I )
232: A( I+1, I ) = ONE
233: *
234: * Compute W(i+1:n,i)
235: *
236: CALL DSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
237: $ A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
238: CALL DGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW,
239: $ A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
240: CALL DGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
241: $ LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
242: CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
243: $ A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
244: CALL DGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
245: $ LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
246: CALL DSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
247: ALPHA = -HALF*TAU( I )*DDOT( N-I, W( I+1, I ), 1,
248: $ A( I+1, I ), 1 )
249: CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
250: END IF
251: *
252: 20 CONTINUE
253: END IF
254: *
255: RETURN
256: *
257: * End of DLATRD
258: *
259: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlatrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack