LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zhptrd.f
Go to the documentation of this file.
1 *> \brief \b ZHPTRD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHPTRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptrd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptrd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptrd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHPTRD( 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 D( * ), E( * )
29 * COMPLEX*16 AP( * ), TAU( * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZHPTRD 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 *> \date November 2011
116 *
117 *> \ingroup complex16OTHERcomputational
118 *
119 *> \par Further Details:
120 * =====================
121 *>
122 *> \verbatim
123 *>
124 *> If UPLO = 'U', the matrix Q is represented as a product of elementary
125 *> reflectors
126 *>
127 *> Q = H(n-1) . . . H(2) H(1).
128 *>
129 *> Each H(i) has the form
130 *>
131 *> H(i) = I - tau * v * v**H
132 *>
133 *> where tau is a complex scalar, and v is a complex vector with
134 *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
135 *> overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
136 *>
137 *> If UPLO = 'L', the matrix Q is represented as a product of elementary
138 *> reflectors
139 *>
140 *> Q = H(1) H(2) . . . H(n-1).
141 *>
142 *> Each H(i) has the form
143 *>
144 *> H(i) = I - tau * v * v**H
145 *>
146 *> where tau is a complex scalar, and v is a complex vector with
147 *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
148 *> overwriting A(i+2:n,i), and tau is stored in TAU(i).
149 *> \endverbatim
150 *>
151 * =====================================================================
152  SUBROUTINE zhptrd( UPLO, N, AP, D, E, TAU, INFO )
153 *
154 * -- LAPACK computational routine (version 3.4.0) --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 * November 2011
158 *
159 * .. Scalar Arguments ..
160  CHARACTER uplo
161  INTEGER info, n
162 * ..
163 * .. Array Arguments ..
164  DOUBLE PRECISION d( * ), e( * )
165  COMPLEX*16 ap( * ), tau( * )
166 * ..
167 *
168 * =====================================================================
169 *
170 * .. Parameters ..
171  COMPLEX*16 one, zero, half
172  parameter( one = ( 1.0d+0, 0.0d+0 ),
173  $ zero = ( 0.0d+0, 0.0d+0 ),
174  $ half = ( 0.5d+0, 0.0d+0 ) )
175 * ..
176 * .. Local Scalars ..
177  LOGICAL upper
178  INTEGER i, i1, i1i1, ii
179  COMPLEX*16 alpha, taui
180 * ..
181 * .. External Subroutines ..
182  EXTERNAL xerbla, zaxpy, zhpmv, zhpr2, zlarfg
183 * ..
184 * .. External Functions ..
185  LOGICAL lsame
186  COMPLEX*16 zdotc
187  EXTERNAL lsame, zdotc
188 * ..
189 * .. Intrinsic Functions ..
190  INTRINSIC dble
191 * ..
192 * .. Executable Statements ..
193 *
194 * Test the input parameters
195 *
196  info = 0
197  upper = lsame( uplo, 'U' )
198  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
199  info = -1
200  ELSE IF( n.LT.0 ) THEN
201  info = -2
202  END IF
203  IF( info.NE.0 ) THEN
204  CALL xerbla( 'ZHPTRD', -info )
205  return
206  END IF
207 *
208 * Quick return if possible
209 *
210  IF( n.LE.0 )
211  $ return
212 *
213  IF( upper ) THEN
214 *
215 * Reduce the upper triangle of A.
216 * I1 is the index in AP of A(1,I+1).
217 *
218  i1 = n*( n-1 ) / 2 + 1
219  ap( i1+n-1 ) = dble( ap( i1+n-1 ) )
220  DO 10 i = n - 1, 1, -1
221 *
222 * Generate elementary reflector H(i) = I - tau * v * v**H
223 * to annihilate A(1:i-1,i+1)
224 *
225  alpha = ap( i1+i-1 )
226  CALL zlarfg( i, alpha, ap( i1 ), 1, taui )
227  e( i ) = alpha
228 *
229  IF( taui.NE.zero ) THEN
230 *
231 * Apply H(i) from both sides to A(1:i,1:i)
232 *
233  ap( i1+i-1 ) = one
234 *
235 * Compute y := tau * A * v storing y in TAU(1:i)
236 *
237  CALL zhpmv( uplo, i, taui, ap, ap( i1 ), 1, zero, tau,
238  $ 1 )
239 *
240 * Compute w := y - 1/2 * tau * (y**H *v) * v
241 *
242  alpha = -half*taui*zdotc( i, tau, 1, ap( i1 ), 1 )
243  CALL zaxpy( i, alpha, ap( i1 ), 1, tau, 1 )
244 *
245 * Apply the transformation as a rank-2 update:
246 * A := A - v * w**H - w * v**H
247 *
248  CALL zhpr2( uplo, i, -one, ap( i1 ), 1, tau, 1, ap )
249 *
250  END IF
251  ap( i1+i-1 ) = e( i )
252  d( i+1 ) = ap( i1+i )
253  tau( i ) = taui
254  i1 = i1 - i
255  10 continue
256  d( 1 ) = ap( 1 )
257  ELSE
258 *
259 * Reduce the lower triangle of A. II is the index in AP of
260 * A(i,i) and I1I1 is the index of A(i+1,i+1).
261 *
262  ii = 1
263  ap( 1 ) = dble( ap( 1 ) )
264  DO 20 i = 1, n - 1
265  i1i1 = ii + n - i + 1
266 *
267 * Generate elementary reflector H(i) = I - tau * v * v**H
268 * to annihilate A(i+2:n,i)
269 *
270  alpha = ap( ii+1 )
271  CALL zlarfg( n-i, alpha, ap( ii+2 ), 1, taui )
272  e( i ) = alpha
273 *
274  IF( taui.NE.zero ) THEN
275 *
276 * Apply H(i) from both sides to A(i+1:n,i+1:n)
277 *
278  ap( ii+1 ) = one
279 *
280 * Compute y := tau * A * v storing y in TAU(i:n-1)
281 *
282  CALL zhpmv( uplo, n-i, taui, ap( i1i1 ), ap( ii+1 ), 1,
283  $ zero, tau( i ), 1 )
284 *
285 * Compute w := y - 1/2 * tau * (y**H *v) * v
286 *
287  alpha = -half*taui*zdotc( n-i, tau( i ), 1, ap( ii+1 ),
288  $ 1 )
289  CALL zaxpy( n-i, alpha, ap( ii+1 ), 1, tau( i ), 1 )
290 *
291 * Apply the transformation as a rank-2 update:
292 * A := A - v * w**H - w * v**H
293 *
294  CALL zhpr2( uplo, n-i, -one, ap( ii+1 ), 1, tau( i ), 1,
295  $ ap( i1i1 ) )
296 *
297  END IF
298  ap( ii+1 ) = e( i )
299  d( i ) = ap( ii )
300  tau( i ) = taui
301  ii = i1i1
302  20 continue
303  d( n ) = ap( ii )
304  END IF
305 *
306  return
307 *
308 * End of ZHPTRD
309 *
310  END