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