|
2329
|
1 SUBROUTINE ZLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) |
|
|
2 * |
|
3333
|
3 * -- LAPACK auxiliary routine (version 3.0) -- |
|
2329
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
6 * June 30, 1999 |
|
2329
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 INTEGER K, LDA, LDT, LDY, N, NB |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ), |
|
|
13 $ Y( LDY, NB ) |
|
|
14 * .. |
|
|
15 * |
|
|
16 * Purpose |
|
|
17 * ======= |
|
|
18 * |
|
|
19 * ZLAHRD reduces the first NB columns of a complex general n-by-(n-k+1) |
|
|
20 * matrix A so that elements below the k-th subdiagonal are zero. The |
|
|
21 * reduction is performed by a unitary similarity transformation |
|
|
22 * Q' * A * Q. The routine returns the matrices V and T which determine |
|
|
23 * Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. |
|
|
24 * |
|
|
25 * This is an auxiliary routine called by ZGEHRD. |
|
|
26 * |
|
|
27 * Arguments |
|
|
28 * ========= |
|
|
29 * |
|
|
30 * N (input) INTEGER |
|
|
31 * The order of the matrix A. |
|
|
32 * |
|
|
33 * K (input) INTEGER |
|
|
34 * The offset for the reduction. Elements below the k-th |
|
|
35 * subdiagonal in the first NB columns are reduced to zero. |
|
|
36 * |
|
|
37 * NB (input) INTEGER |
|
|
38 * The number of columns to be reduced. |
|
|
39 * |
|
|
40 * A (input/output) COMPLEX*16 array, dimension (LDA,N-K+1) |
|
|
41 * On entry, the n-by-(n-k+1) general matrix A. |
|
|
42 * On exit, the elements on and above the k-th subdiagonal in |
|
|
43 * the first NB columns are overwritten with the corresponding |
|
|
44 * elements of the reduced matrix; the elements below the k-th |
|
|
45 * subdiagonal, with the array TAU, represent the matrix Q as a |
|
|
46 * product of elementary reflectors. The other columns of A are |
|
|
47 * unchanged. See Further Details. |
|
|
48 * |
|
|
49 * LDA (input) INTEGER |
|
|
50 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
51 * |
|
|
52 * TAU (output) COMPLEX*16 array, dimension (NB) |
|
|
53 * The scalar factors of the elementary reflectors. See Further |
|
|
54 * Details. |
|
|
55 * |
|
3333
|
56 * T (output) COMPLEX*16 array, dimension (LDT,NB) |
|
2329
|
57 * The upper triangular matrix T. |
|
|
58 * |
|
|
59 * LDT (input) INTEGER |
|
|
60 * The leading dimension of the array T. LDT >= NB. |
|
|
61 * |
|
|
62 * Y (output) COMPLEX*16 array, dimension (LDY,NB) |
|
|
63 * The n-by-nb matrix Y. |
|
|
64 * |
|
|
65 * LDY (input) INTEGER |
|
|
66 * The leading dimension of the array Y. LDY >= max(1,N). |
|
|
67 * |
|
|
68 * Further Details |
|
|
69 * =============== |
|
|
70 * |
|
|
71 * The matrix Q is represented as a product of nb elementary reflectors |
|
|
72 * |
|
|
73 * Q = H(1) H(2) . . . H(nb). |
|
|
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(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in |
|
|
81 * A(i+k+1:n,i), and tau in TAU(i). |
|
|
82 * |
|
|
83 * The elements of the vectors v together form the (n-k+1)-by-nb matrix |
|
|
84 * V which is needed, with T and Y, to apply the transformation to the |
|
|
85 * unreduced part of the matrix, using an update of the form: |
|
|
86 * A := (I - V*T*V') * (A - Y*V'). |
|
|
87 * |
|
|
88 * The contents of A on exit are illustrated by the following example |
|
|
89 * with n = 7, k = 3 and nb = 2: |
|
|
90 * |
|
|
91 * ( a h a a a ) |
|
|
92 * ( a h a a a ) |
|
|
93 * ( a h a a a ) |
|
|
94 * ( h h a a a ) |
|
|
95 * ( v1 h a a a ) |
|
|
96 * ( v1 v2 a a a ) |
|
|
97 * ( v1 v2 a a a ) |
|
|
98 * |
|
|
99 * where a denotes an element of the original matrix A, h denotes a |
|
|
100 * modified element of the upper Hessenberg matrix H, and vi denotes an |
|
|
101 * element of the vector defining H(i). |
|
|
102 * |
|
|
103 * ===================================================================== |
|
|
104 * |
|
|
105 * .. Parameters .. |
|
|
106 COMPLEX*16 ZERO, ONE |
|
|
107 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), |
|
|
108 $ ONE = ( 1.0D+0, 0.0D+0 ) ) |
|
|
109 * .. |
|
|
110 * .. Local Scalars .. |
|
|
111 INTEGER I |
|
|
112 COMPLEX*16 EI |
|
|
113 * .. |
|
|
114 * .. External Subroutines .. |
|
|
115 EXTERNAL ZAXPY, ZCOPY, ZGEMV, ZLACGV, ZLARFG, ZSCAL, |
|
|
116 $ ZTRMV |
|
|
117 * .. |
|
|
118 * .. Intrinsic Functions .. |
|
|
119 INTRINSIC MIN |
|
|
120 * .. |
|
|
121 * .. Executable Statements .. |
|
|
122 * |
|
|
123 * Quick return if possible |
|
|
124 * |
|
|
125 IF( N.LE.1 ) |
|
|
126 $ RETURN |
|
|
127 * |
|
|
128 DO 10 I = 1, NB |
|
|
129 IF( I.GT.1 ) THEN |
|
|
130 * |
|
|
131 * Update A(1:n,i) |
|
|
132 * |
|
|
133 * Compute i-th column of A - Y * V' |
|
|
134 * |
|
|
135 CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) |
|
|
136 CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, |
|
|
137 $ A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 ) |
|
|
138 CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA ) |
|
|
139 * |
|
|
140 * Apply I - V * T' * V' to this column (call it b) from the |
|
|
141 * left, using the last column of T as workspace |
|
|
142 * |
|
|
143 * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) |
|
|
144 * ( V2 ) ( b2 ) |
|
|
145 * |
|
|
146 * where V1 is unit lower triangular |
|
|
147 * |
|
|
148 * w := V1' * b1 |
|
|
149 * |
|
|
150 CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) |
|
|
151 CALL ZTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1, |
|
|
152 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) |
|
|
153 * |
|
|
154 * w := w + V2'*b2 |
|
|
155 * |
|
|
156 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE, |
|
|
157 $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE, |
|
|
158 $ T( 1, NB ), 1 ) |
|
|
159 * |
|
|
160 * w := T'*w |
|
|
161 * |
|
|
162 CALL ZTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1, |
|
|
163 $ T, LDT, T( 1, NB ), 1 ) |
|
|
164 * |
|
|
165 * b2 := b2 - V2*w |
|
|
166 * |
|
|
167 CALL ZGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ), |
|
|
168 $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) |
|
|
169 * |
|
|
170 * b1 := b1 - V1*w |
|
|
171 * |
|
|
172 CALL ZTRMV( 'Lower', 'No transpose', 'Unit', I-1, |
|
|
173 $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) |
|
|
174 CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) |
|
|
175 * |
|
|
176 A( K+I-1, I-1 ) = EI |
|
|
177 END IF |
|
|
178 * |
|
|
179 * Generate the elementary reflector H(i) to annihilate |
|
|
180 * A(k+i+1:n,i) |
|
|
181 * |
|
|
182 EI = A( K+I, I ) |
|
|
183 CALL ZLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1, |
|
|
184 $ TAU( I ) ) |
|
|
185 A( K+I, I ) = ONE |
|
|
186 * |
|
|
187 * Compute Y(1:n,i) |
|
|
188 * |
|
|
189 CALL ZGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA, |
|
|
190 $ A( K+I, I ), 1, ZERO, Y( 1, I ), 1 ) |
|
|
191 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE, |
|
|
192 $ A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ), |
|
|
193 $ 1 ) |
|
|
194 CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1, |
|
|
195 $ ONE, Y( 1, I ), 1 ) |
|
|
196 CALL ZSCAL( N, TAU( I ), Y( 1, I ), 1 ) |
|
|
197 * |
|
|
198 * Compute T(1:i,i) |
|
|
199 * |
|
|
200 CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) |
|
|
201 CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT, |
|
|
202 $ T( 1, I ), 1 ) |
|
|
203 T( I, I ) = TAU( I ) |
|
|
204 * |
|
|
205 10 CONTINUE |
|
|
206 A( K+NB, NB ) = EI |
|
|
207 * |
|
|
208 RETURN |
|
|
209 * |
|
|
210 * End of ZLAHRD |
|
|
211 * |
|
|
212 END |
