LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsptrd.f
Go to the documentation of this file.
1*> \brief \b DSPTRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSPTRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSPTRD( UPLO, N, AP, D, E, TAU, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, N
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION AP( * ), D( * ), E( * ), TAU( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DSPTRD reduces a real symmetric matrix A stored in packed form to
38*> symmetric tridiagonal form T by an orthogonal similarity
39*> transformation: Q**T * A * Q = T.
40*> \endverbatim
41*
42* Arguments:
43* ==========
44*
45*> \param[in] UPLO
46*> \verbatim
47*> UPLO is CHARACTER*1
48*> = 'U': Upper triangle of A is stored;
49*> = 'L': Lower triangle of A is stored.
50*> \endverbatim
51*>
52*> \param[in] N
53*> \verbatim
54*> N is INTEGER
55*> The order of the matrix A. N >= 0.
56*> \endverbatim
57*>
58*> \param[in,out] AP
59*> \verbatim
60*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
61*> On entry, the upper or lower triangle of the symmetric matrix
62*> A, packed columnwise in a linear array. The j-th column of A
63*> is stored in the array AP as follows:
64*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
65*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
66*> On exit, if UPLO = 'U', the diagonal and first superdiagonal
67*> of A are overwritten by the corresponding elements of the
68*> tridiagonal matrix T, and the elements above the first
69*> superdiagonal, with the array TAU, represent the orthogonal
70*> matrix Q as a product of elementary reflectors; if UPLO
71*> = 'L', the diagonal and first subdiagonal of A are over-
72*> written by the corresponding elements of the tridiagonal
73*> matrix T, and the elements below the first subdiagonal, with
74*> the array TAU, represent the orthogonal matrix Q as a product
75*> of elementary reflectors. See Further Details.
76*> \endverbatim
77*>
78*> \param[out] D
79*> \verbatim
80*> D is DOUBLE PRECISION array, dimension (N)
81*> The diagonal elements of the tridiagonal matrix T:
82*> D(i) = A(i,i).
83*> \endverbatim
84*>
85*> \param[out] E
86*> \verbatim
87*> E is DOUBLE PRECISION array, dimension (N-1)
88*> The off-diagonal elements of the tridiagonal matrix T:
89*> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
90*> \endverbatim
91*>
92*> \param[out] TAU
93*> \verbatim
94*> TAU is DOUBLE PRECISION array, dimension (N-1)
95*> The scalar factors of the elementary reflectors (see Further
96*> Details).
97*> \endverbatim
98*>
99*> \param[out] INFO
100*> \verbatim
101*> INFO is INTEGER
102*> = 0: successful exit
103*> < 0: if INFO = -i, the i-th argument had an illegal value
104*> \endverbatim
105*
106* Authors:
107* ========
108*
109*> \author Univ. of Tennessee
110*> \author Univ. of California Berkeley
111*> \author Univ. of Colorado Denver
112*> \author NAG Ltd.
113*
114*> \ingroup hptrd
115*
116*> \par Further Details:
117* =====================
118*>
119*> \verbatim
120*>
121*> If UPLO = 'U', the matrix Q is represented as a product of elementary
122*> reflectors
123*>
124*> Q = H(n-1) . . . H(2) H(1).
125*>
126*> Each H(i) has the form
127*>
128*> H(i) = I - tau * v * v**T
129*>
130*> where tau is a real scalar, and v is a real vector with
131*> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
132*> overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
133*>
134*> If UPLO = 'L', the matrix Q is represented as a product of elementary
135*> reflectors
136*>
137*> Q = H(1) H(2) . . . H(n-1).
138*>
139*> Each H(i) has the form
140*>
141*> H(i) = I - tau * v * v**T
142*>
143*> where tau is a real scalar, and v is a real vector with
144*> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
145*> overwriting A(i+2:n,i), and tau is stored in TAU(i).
146*> \endverbatim
147*>
148* =====================================================================
149 SUBROUTINE dsptrd( UPLO, N, AP, D, E, TAU, INFO )
150*
151* -- LAPACK computational routine --
152* -- LAPACK is a software package provided by Univ. of Tennessee, --
153* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
154*
155* .. Scalar Arguments ..
156 CHARACTER UPLO
157 INTEGER INFO, N
158* ..
159* .. Array Arguments ..
160 DOUBLE PRECISION AP( * ), D( * ), E( * ), TAU( * )
161* ..
162*
163* =====================================================================
164*
165* .. Parameters ..
166 DOUBLE PRECISION ONE, ZERO, HALF
167 parameter( one = 1.0d0, zero = 0.0d0,
168 $ half = 1.0d0 / 2.0d0 )
169* ..
170* .. Local Scalars ..
171 LOGICAL UPPER
172 INTEGER I, I1, I1I1, II
173 DOUBLE PRECISION ALPHA, TAUI
174* ..
175* .. External Subroutines ..
176 EXTERNAL daxpy, dlarfg, dspmv, dspr2, xerbla
177* ..
178* .. External Functions ..
179 LOGICAL LSAME
180 DOUBLE PRECISION DDOT
181 EXTERNAL lsame, ddot
182* ..
183* .. Executable Statements ..
184*
185* Test the input parameters
186*
187 info = 0
188 upper = lsame( uplo, 'U' )
189 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
190 info = -1
191 ELSE IF( n.LT.0 ) THEN
192 info = -2
193 END IF
194 IF( info.NE.0 ) THEN
195 CALL xerbla( 'DSPTRD', -info )
196 RETURN
197 END IF
198*
199* Quick return if possible
200*
201 IF( n.LE.0 )
202 $ RETURN
203*
204 IF( upper ) THEN
205*
206* Reduce the upper triangle of A.
207* I1 is the index in AP of A(1,I+1).
208*
209 i1 = n*( n-1 ) / 2 + 1
210 DO 10 i = n - 1, 1, -1
211*
212* Generate elementary reflector H(i) = I - tau * v * v**T
213* to annihilate A(1:i-1,i+1)
214*
215 CALL dlarfg( i, ap( i1+i-1 ), ap( i1 ), 1, taui )
216 e( i ) = ap( i1+i-1 )
217*
218 IF( taui.NE.zero ) THEN
219*
220* Apply H(i) from both sides to A(1:i,1:i)
221*
222 ap( i1+i-1 ) = one
223*
224* Compute y := tau * A * v storing y in TAU(1:i)
225*
226 CALL dspmv( uplo, i, taui, ap, ap( i1 ), 1, zero, tau,
227 $ 1 )
228*
229* Compute w := y - 1/2 * tau * (y**T *v) * v
230*
231 alpha = -half*taui*ddot( i, tau, 1, ap( i1 ), 1 )
232 CALL daxpy( i, alpha, ap( i1 ), 1, tau, 1 )
233*
234* Apply the transformation as a rank-2 update:
235* A := A - v * w**T - w * v**T
236*
237 CALL dspr2( uplo, i, -one, ap( i1 ), 1, tau, 1, ap )
238*
239 ap( i1+i-1 ) = e( i )
240 END IF
241 d( i+1 ) = ap( i1+i )
242 tau( i ) = taui
243 i1 = i1 - i
244 10 CONTINUE
245 d( 1 ) = ap( 1 )
246 ELSE
247*
248* Reduce the lower triangle of A. II is the index in AP of
249* A(i,i) and I1I1 is the index of A(i+1,i+1).
250*
251 ii = 1
252 DO 20 i = 1, n - 1
253 i1i1 = ii + n - i + 1
254*
255* Generate elementary reflector H(i) = I - tau * v * v**T
256* to annihilate A(i+2:n,i)
257*
258 CALL dlarfg( n-i, ap( ii+1 ), ap( ii+2 ), 1, taui )
259 e( i ) = ap( ii+1 )
260*
261 IF( taui.NE.zero ) THEN
262*
263* Apply H(i) from both sides to A(i+1:n,i+1:n)
264*
265 ap( ii+1 ) = one
266*
267* Compute y := tau * A * v storing y in TAU(i:n-1)
268*
269 CALL dspmv( uplo, n-i, taui, ap( i1i1 ), ap( ii+1 ), 1,
270 $ zero, tau( i ), 1 )
271*
272* Compute w := y - 1/2 * tau * (y**T *v) * v
273*
274 alpha = -half*taui*ddot( n-i, tau( i ), 1, ap( ii+1 ),
275 $ 1 )
276 CALL daxpy( n-i, alpha, ap( ii+1 ), 1, tau( i ), 1 )
277*
278* Apply the transformation as a rank-2 update:
279* A := A - v * w**T - w * v**T
280*
281 CALL dspr2( uplo, n-i, -one, ap( ii+1 ), 1, tau( i ), 1,
282 $ ap( i1i1 ) )
283*
284 ap( ii+1 ) = e( i )
285 END IF
286 d( i ) = ap( ii )
287 tau( i ) = taui
288 ii = i1i1
289 20 CONTINUE
290 d( n ) = ap( ii )
291 END IF
292*
293 RETURN
294*
295* End of DSPTRD
296*
297 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dspmv(uplo, n, alpha, ap, x, incx, beta, y, incy)
DSPMV
Definition dspmv.f:147
subroutine dspr2(uplo, n, alpha, x, incx, y, incy, ap)
DSPR2
Definition dspr2.f:142
subroutine dsptrd(uplo, n, ap, d, e, tau, info)
DSPTRD
Definition dsptrd.f:150
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106