Annotation of rpl/lapack/lapack/zhptrd.f, revision 1.3
1.1 bertrand 1: SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER INFO, N
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION D( * ), E( * )
14: COMPLEX*16 AP( * ), TAU( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZHPTRD reduces a complex Hermitian matrix A stored in packed form to
21: * real symmetric tridiagonal form T by a unitary similarity
22: * transformation: Q**H * A * Q = T.
23: *
24: * Arguments
25: * =========
26: *
27: * UPLO (input) CHARACTER*1
28: * = 'U': Upper triangle of A is stored;
29: * = 'L': Lower triangle of A is stored.
30: *
31: * N (input) INTEGER
32: * The order of the matrix A. N >= 0.
33: *
34: * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
35: * On entry, the upper or lower triangle of the Hermitian matrix
36: * A, packed columnwise in a linear array. The j-th column of A
37: * is stored in the array AP as follows:
38: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
39: * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
40: * On exit, if UPLO = 'U', the diagonal and first superdiagonal
41: * of A are overwritten by the corresponding elements of the
42: * tridiagonal matrix T, and the elements above the first
43: * superdiagonal, with the array TAU, represent the unitary
44: * matrix Q as a product of elementary reflectors; if UPLO
45: * = 'L', the diagonal and first subdiagonal of A are over-
46: * written by the corresponding elements of the tridiagonal
47: * matrix T, and the elements below the first subdiagonal, with
48: * the array TAU, represent the unitary matrix Q as a product
49: * of elementary reflectors. See Further Details.
50: *
51: * D (output) DOUBLE PRECISION array, dimension (N)
52: * The diagonal elements of the tridiagonal matrix T:
53: * D(i) = A(i,i).
54: *
55: * E (output) DOUBLE PRECISION array, dimension (N-1)
56: * The off-diagonal elements of the tridiagonal matrix T:
57: * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
58: *
59: * TAU (output) COMPLEX*16 array, dimension (N-1)
60: * The scalar factors of the elementary reflectors (see Further
61: * Details).
62: *
63: * INFO (output) INTEGER
64: * = 0: successful exit
65: * < 0: if INFO = -i, the i-th argument had an illegal value
66: *
67: * Further Details
68: * ===============
69: *
70: * If UPLO = 'U', the matrix Q is represented as a product of elementary
71: * reflectors
72: *
73: * Q = H(n-1) . . . H(2) H(1).
74: *
75: * Each H(i) has the form
76: *
77: * H(i) = I - tau * v * v'
78: *
79: * where tau is a complex scalar, and v is a complex vector with
80: * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
81: * overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
82: *
83: * If UPLO = 'L', the matrix Q is represented as a product of elementary
84: * reflectors
85: *
86: * Q = H(1) H(2) . . . H(n-1).
87: *
88: * Each H(i) has the form
89: *
90: * H(i) = I - tau * v * v'
91: *
92: * where tau is a complex scalar, and v is a complex vector with
93: * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
94: * overwriting A(i+2:n,i), and tau is stored in TAU(i).
95: *
96: * =====================================================================
97: *
98: * .. Parameters ..
99: COMPLEX*16 ONE, ZERO, HALF
100: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
101: $ ZERO = ( 0.0D+0, 0.0D+0 ),
102: $ HALF = ( 0.5D+0, 0.0D+0 ) )
103: * ..
104: * .. Local Scalars ..
105: LOGICAL UPPER
106: INTEGER I, I1, I1I1, II
107: COMPLEX*16 ALPHA, TAUI
108: * ..
109: * .. External Subroutines ..
110: EXTERNAL XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG
111: * ..
112: * .. External Functions ..
113: LOGICAL LSAME
114: COMPLEX*16 ZDOTC
115: EXTERNAL LSAME, ZDOTC
116: * ..
117: * .. Intrinsic Functions ..
118: INTRINSIC DBLE
119: * ..
120: * .. Executable Statements ..
121: *
122: * Test the input parameters
123: *
124: INFO = 0
125: UPPER = LSAME( UPLO, 'U' )
126: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
127: INFO = -1
128: ELSE IF( N.LT.0 ) THEN
129: INFO = -2
130: END IF
131: IF( INFO.NE.0 ) THEN
132: CALL XERBLA( 'ZHPTRD', -INFO )
133: RETURN
134: END IF
135: *
136: * Quick return if possible
137: *
138: IF( N.LE.0 )
139: $ RETURN
140: *
141: IF( UPPER ) THEN
142: *
143: * Reduce the upper triangle of A.
144: * I1 is the index in AP of A(1,I+1).
145: *
146: I1 = N*( N-1 ) / 2 + 1
147: AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) )
148: DO 10 I = N - 1, 1, -1
149: *
150: * Generate elementary reflector H(i) = I - tau * v * v'
151: * to annihilate A(1:i-1,i+1)
152: *
153: ALPHA = AP( I1+I-1 )
154: CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI )
155: E( I ) = ALPHA
156: *
157: IF( TAUI.NE.ZERO ) THEN
158: *
159: * Apply H(i) from both sides to A(1:i,1:i)
160: *
161: AP( I1+I-1 ) = ONE
162: *
163: * Compute y := tau * A * v storing y in TAU(1:i)
164: *
165: CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
166: $ 1 )
167: *
168: * Compute w := y - 1/2 * tau * (y'*v) * v
169: *
170: ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 )
171: CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
172: *
173: * Apply the transformation as a rank-2 update:
174: * A := A - v * w' - w * v'
175: *
176: CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
177: *
178: END IF
179: AP( I1+I-1 ) = E( I )
180: D( I+1 ) = AP( I1+I )
181: TAU( I ) = TAUI
182: I1 = I1 - I
183: 10 CONTINUE
184: D( 1 ) = AP( 1 )
185: ELSE
186: *
187: * Reduce the lower triangle of A. II is the index in AP of
188: * A(i,i) and I1I1 is the index of A(i+1,i+1).
189: *
190: II = 1
191: AP( 1 ) = DBLE( AP( 1 ) )
192: DO 20 I = 1, N - 1
193: I1I1 = II + N - I + 1
194: *
195: * Generate elementary reflector H(i) = I - tau * v * v'
196: * to annihilate A(i+2:n,i)
197: *
198: ALPHA = AP( II+1 )
199: CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI )
200: E( I ) = ALPHA
201: *
202: IF( TAUI.NE.ZERO ) THEN
203: *
204: * Apply H(i) from both sides to A(i+1:n,i+1:n)
205: *
206: AP( II+1 ) = ONE
207: *
208: * Compute y := tau * A * v storing y in TAU(i:n-1)
209: *
210: CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
211: $ ZERO, TAU( I ), 1 )
212: *
213: * Compute w := y - 1/2 * tau * (y'*v) * v
214: *
215: ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ),
216: $ 1 )
217: CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
218: *
219: * Apply the transformation as a rank-2 update:
220: * A := A - v * w' - w * v'
221: *
222: CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
223: $ AP( I1I1 ) )
224: *
225: END IF
226: AP( II+1 ) = E( I )
227: D( I ) = AP( II )
228: TAU( I ) = TAUI
229: II = I1I1
230: 20 CONTINUE
231: D( N ) = AP( II )
232: END IF
233: *
234: RETURN
235: *
236: * End of ZHPTRD
237: *
238: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhptrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack