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*> \date November 2011
116*
117*> \ingroup complexOTHERcomputational
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 CHPTRD( 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      REAL               D( * ), E( * )
165      COMPLEX            AP( * ), TAU( * )
166*     ..
167*
168*  =====================================================================
169*
170*     .. Parameters ..
171      COMPLEX            ONE, ZERO, HALF
172      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
173     $                   ZERO = ( 0.0E+0, 0.0E+0 ),
174     $                   HALF = ( 0.5E+0, 0.0E+0 ) )
175*     ..
176*     .. Local Scalars ..
177      LOGICAL            UPPER
178      INTEGER            I, I1, I1I1, II
179      COMPLEX            ALPHA, TAUI
180*     ..
181*     .. External Subroutines ..
182      EXTERNAL           CAXPY, CHPMV, CHPR2, CLARFG, XERBLA
183*     ..
184*     .. External Functions ..
185      LOGICAL            LSAME
186      COMPLEX            CDOTC
187      EXTERNAL           LSAME, CDOTC
188*     ..
189*     .. Intrinsic Functions ..
190      INTRINSIC          REAL
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( 'CHPTRD', -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 ) = REAL( 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 CLARFG( 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 CHPMV( 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*CDOTC( I, TAU, 1, AP( I1 ), 1 )
243               CALL CAXPY( 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 CHPR2( 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 ) = REAL( 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 CLARFG( 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 CHPMV( 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*CDOTC( N-I, TAU( I ), 1, AP( II+1 ),
288     $                 1 )
289               CALL CAXPY( 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 CHPR2( 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 CHPTRD
309*
310      END
311