1 /* ./src_f77/clahrd.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static complex c_b1 = {0.f,0.f};
11 static complex c_b2 = {1.f,0.f};
12 static integer c__1 = 1;
13
clahrd_(integer * n,integer * k,integer * nb,complex * a,integer * lda,complex * tau,complex * t,integer * ldt,complex * y,integer * ldy)14 /* Subroutine */ int clahrd_(integer *n, integer *k, integer *nb, complex *a,
15 integer *lda, complex *tau, complex *t, integer *ldt, complex *y,
16 integer *ldy)
17 {
18 /* System generated locals */
19 integer a_dim1, a_offset, t_dim1, t_offset, y_dim1, y_offset, i__1, i__2,
20 i__3;
21 complex q__1;
22
23 /* Local variables */
24 static integer i__;
25 static complex ei;
26 extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
27 integer *), cgemv_(char *, integer *, integer *, complex *,
28 complex *, integer *, complex *, integer *, complex *, complex *,
29 integer *, ftnlen), ccopy_(integer *, complex *, integer *,
30 complex *, integer *), caxpy_(integer *, complex *, complex *,
31 integer *, complex *, integer *), ctrmv_(char *, char *, char *,
32 integer *, complex *, integer *, complex *, integer *, ftnlen,
33 ftnlen, ftnlen), clarfg_(integer *, complex *, complex *, integer
34 *, complex *), clacgv_(integer *, complex *, integer *);
35
36
37 /* -- LAPACK auxiliary routine (version 3.0) -- */
38 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
39 /* Courant Institute, Argonne National Lab, and Rice University */
40 /* June 30, 1999 */
41
42 /* .. Scalar Arguments .. */
43 /* .. */
44 /* .. Array Arguments .. */
45 /* .. */
46
47 /* Purpose */
48 /* ======= */
49
50 /* CLAHRD reduces the first NB columns of a complex general n-by-(n-k+1) */
51 /* matrix A so that elements below the k-th subdiagonal are zero. The */
52 /* reduction is performed by a unitary similarity transformation */
53 /* Q' * A * Q. The routine returns the matrices V and T which determine */
54 /* Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. */
55
56 /* This is an auxiliary routine called by CGEHRD. */
57
58 /* Arguments */
59 /* ========= */
60
61 /* N (input) INTEGER */
62 /* The order of the matrix A. */
63
64 /* K (input) INTEGER */
65 /* The offset for the reduction. Elements below the k-th */
66 /* subdiagonal in the first NB columns are reduced to zero. */
67
68 /* NB (input) INTEGER */
69 /* The number of columns to be reduced. */
70
71 /* A (input/output) COMPLEX array, dimension (LDA,N-K+1) */
72 /* On entry, the n-by-(n-k+1) general matrix A. */
73 /* On exit, the elements on and above the k-th subdiagonal in */
74 /* the first NB columns are overwritten with the corresponding */
75 /* elements of the reduced matrix; the elements below the k-th */
76 /* subdiagonal, with the array TAU, represent the matrix Q as a */
77 /* product of elementary reflectors. The other columns of A are */
78 /* unchanged. See Further Details. */
79
80 /* LDA (input) INTEGER */
81 /* The leading dimension of the array A. LDA >= max(1,N). */
82
83 /* TAU (output) COMPLEX array, dimension (NB) */
84 /* The scalar factors of the elementary reflectors. See Further */
85 /* Details. */
86
87 /* T (output) COMPLEX array, dimension (LDT,NB) */
88 /* The upper triangular matrix T. */
89
90 /* LDT (input) INTEGER */
91 /* The leading dimension of the array T. LDT >= NB. */
92
93 /* Y (output) COMPLEX array, dimension (LDY,NB) */
94 /* The n-by-nb matrix Y. */
95
96 /* LDY (input) INTEGER */
97 /* The leading dimension of the array Y. LDY >= max(1,N). */
98
99 /* Further Details */
100 /* =============== */
101
102 /* The matrix Q is represented as a product of nb elementary reflectors */
103
104 /* Q = H(1) H(2) . . . H(nb). */
105
106 /* Each H(i) has the form */
107
108 /* H(i) = I - tau * v * v' */
109
110 /* where tau is a complex scalar, and v is a complex vector with */
111 /* v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in */
112 /* A(i+k+1:n,i), and tau in TAU(i). */
113
114 /* The elements of the vectors v together form the (n-k+1)-by-nb matrix */
115 /* V which is needed, with T and Y, to apply the transformation to the */
116 /* unreduced part of the matrix, using an update of the form: */
117 /* A := (I - V*T*V') * (A - Y*V'). */
118
119 /* The contents of A on exit are illustrated by the following example */
120 /* with n = 7, k = 3 and nb = 2: */
121
122 /* ( a h a a a ) */
123 /* ( a h a a a ) */
124 /* ( a h a a a ) */
125 /* ( h h a a a ) */
126 /* ( v1 h a a a ) */
127 /* ( v1 v2 a a a ) */
128 /* ( v1 v2 a a a ) */
129
130 /* where a denotes an element of the original matrix A, h denotes a */
131 /* modified element of the upper Hessenberg matrix H, and vi denotes an */
132 /* element of the vector defining H(i). */
133
134 /* ===================================================================== */
135
136 /* .. Parameters .. */
137 /* .. */
138 /* .. Local Scalars .. */
139 /* .. */
140 /* .. External Subroutines .. */
141 /* .. */
142 /* .. Intrinsic Functions .. */
143 /* .. */
144 /* .. Executable Statements .. */
145
146 /* Quick return if possible */
147
148 /* Parameter adjustments */
149 --tau;
150 a_dim1 = *lda;
151 a_offset = 1 + a_dim1;
152 a -= a_offset;
153 t_dim1 = *ldt;
154 t_offset = 1 + t_dim1;
155 t -= t_offset;
156 y_dim1 = *ldy;
157 y_offset = 1 + y_dim1;
158 y -= y_offset;
159
160 /* Function Body */
161 if (*n <= 1) {
162 return 0;
163 }
164
165 i__1 = *nb;
166 for (i__ = 1; i__ <= i__1; ++i__) {
167 if (i__ > 1) {
168
169 /* Update A(1:n,i) */
170
171 /* Compute i-th column of A - Y * V' */
172
173 i__2 = i__ - 1;
174 clacgv_(&i__2, &a[*k + i__ - 1 + a_dim1], lda);
175 i__2 = i__ - 1;
176 q__1.r = -1.f, q__1.i = -0.f;
177 cgemv_("No transpose", n, &i__2, &q__1, &y[y_offset], ldy, &a[*k
178 + i__ - 1 + a_dim1], lda, &c_b2, &a[i__ * a_dim1 + 1], &
179 c__1, (ftnlen)12);
180 i__2 = i__ - 1;
181 clacgv_(&i__2, &a[*k + i__ - 1 + a_dim1], lda);
182
183 /* Apply I - V * T' * V' to this column (call it b) from the */
184 /* left, using the last column of T as workspace */
185
186 /* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) */
187 /* ( V2 ) ( b2 ) */
188
189 /* where V1 is unit lower triangular */
190
191 /* w := V1' * b1 */
192
193 i__2 = i__ - 1;
194 ccopy_(&i__2, &a[*k + 1 + i__ * a_dim1], &c__1, &t[*nb * t_dim1 +
195 1], &c__1);
196 i__2 = i__ - 1;
197 ctrmv_("Lower", "Conjugate transpose", "Unit", &i__2, &a[*k + 1 +
198 a_dim1], lda, &t[*nb * t_dim1 + 1], &c__1, (ftnlen)5, (
199 ftnlen)19, (ftnlen)4);
200
201 /* w := w + V2'*b2 */
202
203 i__2 = *n - *k - i__ + 1;
204 i__3 = i__ - 1;
205 cgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i__ +
206 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b2, &
207 t[*nb * t_dim1 + 1], &c__1, (ftnlen)19);
208
209 /* w := T'*w */
210
211 i__2 = i__ - 1;
212 ctrmv_("Upper", "Conjugate transpose", "Non-unit", &i__2, &t[
213 t_offset], ldt, &t[*nb * t_dim1 + 1], &c__1, (ftnlen)5, (
214 ftnlen)19, (ftnlen)8);
215
216 /* b2 := b2 - V2*w */
217
218 i__2 = *n - *k - i__ + 1;
219 i__3 = i__ - 1;
220 q__1.r = -1.f, q__1.i = -0.f;
221 cgemv_("No transpose", &i__2, &i__3, &q__1, &a[*k + i__ + a_dim1],
222 lda, &t[*nb * t_dim1 + 1], &c__1, &c_b2, &a[*k + i__ +
223 i__ * a_dim1], &c__1, (ftnlen)12);
224
225 /* b1 := b1 - V1*w */
226
227 i__2 = i__ - 1;
228 ctrmv_("Lower", "No transpose", "Unit", &i__2, &a[*k + 1 + a_dim1]
229 , lda, &t[*nb * t_dim1 + 1], &c__1, (ftnlen)5, (ftnlen)12,
230 (ftnlen)4);
231 i__2 = i__ - 1;
232 q__1.r = -1.f, q__1.i = -0.f;
233 caxpy_(&i__2, &q__1, &t[*nb * t_dim1 + 1], &c__1, &a[*k + 1 + i__
234 * a_dim1], &c__1);
235
236 i__2 = *k + i__ - 1 + (i__ - 1) * a_dim1;
237 a[i__2].r = ei.r, a[i__2].i = ei.i;
238 }
239
240 /* Generate the elementary reflector H(i) to annihilate */
241 /* A(k+i+1:n,i) */
242
243 i__2 = *k + i__ + i__ * a_dim1;
244 ei.r = a[i__2].r, ei.i = a[i__2].i;
245 i__2 = *n - *k - i__ + 1;
246 /* Computing MIN */
247 i__3 = *k + i__ + 1;
248 clarfg_(&i__2, &ei, &a[min(i__3,*n) + i__ * a_dim1], &c__1, &tau[i__])
249 ;
250 i__2 = *k + i__ + i__ * a_dim1;
251 a[i__2].r = 1.f, a[i__2].i = 0.f;
252
253 /* Compute Y(1:n,i) */
254
255 i__2 = *n - *k - i__ + 1;
256 cgemv_("No transpose", n, &i__2, &c_b2, &a[(i__ + 1) * a_dim1 + 1],
257 lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b1, &y[i__ *
258 y_dim1 + 1], &c__1, (ftnlen)12);
259 i__2 = *n - *k - i__ + 1;
260 i__3 = i__ - 1;
261 cgemv_("Conjugate transpose", &i__2, &i__3, &c_b2, &a[*k + i__ +
262 a_dim1], lda, &a[*k + i__ + i__ * a_dim1], &c__1, &c_b1, &t[
263 i__ * t_dim1 + 1], &c__1, (ftnlen)19);
264 i__2 = i__ - 1;
265 q__1.r = -1.f, q__1.i = -0.f;
266 cgemv_("No transpose", n, &i__2, &q__1, &y[y_offset], ldy, &t[i__ *
267 t_dim1 + 1], &c__1, &c_b2, &y[i__ * y_dim1 + 1], &c__1, (
268 ftnlen)12);
269 cscal_(n, &tau[i__], &y[i__ * y_dim1 + 1], &c__1);
270
271 /* Compute T(1:i,i) */
272
273 i__2 = i__ - 1;
274 i__3 = i__;
275 q__1.r = -tau[i__3].r, q__1.i = -tau[i__3].i;
276 cscal_(&i__2, &q__1, &t[i__ * t_dim1 + 1], &c__1);
277 i__2 = i__ - 1;
278 ctrmv_("Upper", "No transpose", "Non-unit", &i__2, &t[t_offset], ldt,
279 &t[i__ * t_dim1 + 1], &c__1, (ftnlen)5, (ftnlen)12, (ftnlen)8)
280 ;
281 i__2 = i__ + i__ * t_dim1;
282 i__3 = i__;
283 t[i__2].r = tau[i__3].r, t[i__2].i = tau[i__3].i;
284
285 /* L10: */
286 }
287 i__1 = *k + *nb + *nb * a_dim1;
288 a[i__1].r = ei.r, a[i__1].i = ei.i;
289
290 return 0;
291
292 /* End of CLAHRD */
293
294 } /* clahrd_ */
295
296