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