LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 doubleOTHERcomputational
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)
XERBLA
Definition: xerbla.f:60
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
DSPR2
Definition: dspr2.f:142
subroutine dspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
DSPMV
Definition: dspmv.f:147
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106
subroutine dsptrd(UPLO, N, AP, D, E, TAU, INFO)
DSPTRD
Definition: dsptrd.f:150