1 /* ./src_f77/cunmrz.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 integer c__1 = 1;
11 static integer c_n1 = -1;
12 static integer c__2 = 2;
13 static integer c__65 = 65;
14
cunmrz_(char * side,char * trans,integer * m,integer * n,integer * k,integer * l,complex * a,integer * lda,complex * tau,complex * c__,integer * ldc,complex * work,integer * lwork,integer * info,ftnlen side_len,ftnlen trans_len)15 /* Subroutine */ int cunmrz_(char *side, char *trans, integer *m, integer *n,
16 integer *k, integer *l, complex *a, integer *lda, complex *tau,
17 complex *c__, integer *ldc, complex *work, integer *lwork, integer *
18 info, ftnlen side_len, ftnlen trans_len)
19 {
20 /* System generated locals */
21 address a__1[2];
22 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3[2], i__4,
23 i__5;
24 char ch__1[2];
25
26 /* Builtin functions */
27 /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
28
29 /* Local variables */
30 static integer i__;
31 static complex t[4160] /* was [65][64] */;
32 static integer i1, i2, i3, ib, ic, ja, jc, nb, mi, ni, nq, nw, iws;
33 static logical left;
34 extern logical lsame_(char *, char *, ftnlen, ftnlen);
35 static integer nbmin, iinfo;
36 extern /* Subroutine */ int cunmr3_(char *, char *, integer *, integer *,
37 integer *, integer *, complex *, integer *, complex *, complex *,
38 integer *, complex *, integer *, ftnlen, ftnlen), clarzb_(char *,
39 char *, char *, char *, integer *, integer *, integer *, integer *
40 , complex *, integer *, complex *, integer *, complex *, integer *
41 , complex *, integer *, ftnlen, ftnlen, ftnlen, ftnlen);
42 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
43 integer *, integer *, ftnlen, ftnlen);
44 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), clarzt_(
45 char *, char *, integer *, integer *, complex *, integer *,
46 complex *, complex *, integer *, ftnlen, ftnlen);
47 static logical notran;
48 static integer ldwork;
49 static char transt[1];
50 static integer lwkopt;
51 static logical lquery;
52
53
54 /* -- LAPACK routine (version 3.0) -- */
55 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
56 /* Courant Institute, Argonne National Lab, and Rice University */
57 /* June 30, 1999 */
58
59 /* .. Scalar Arguments .. */
60 /* .. */
61 /* .. Array Arguments .. */
62 /* .. */
63
64 /* Purpose */
65 /* ======= */
66
67 /* CUNMRZ overwrites the general complex M-by-N matrix C with */
68
69 /* SIDE = 'L' SIDE = 'R' */
70 /* TRANS = 'N': Q * C C * Q */
71 /* TRANS = 'C': Q**H * C C * Q**H */
72
73 /* where Q is a complex unitary matrix defined as the product of k */
74 /* elementary reflectors */
75
76 /* Q = H(1) H(2) . . . H(k) */
77
78 /* as returned by CTZRZF. Q is of order M if SIDE = 'L' and of order N */
79 /* if SIDE = 'R'. */
80
81 /* Arguments */
82 /* ========= */
83
84 /* SIDE (input) CHARACTER*1 */
85 /* = 'L': apply Q or Q**H from the Left; */
86 /* = 'R': apply Q or Q**H from the Right. */
87
88 /* TRANS (input) CHARACTER*1 */
89 /* = 'N': No transpose, apply Q; */
90 /* = 'C': Conjugate transpose, apply Q**H. */
91
92 /* M (input) INTEGER */
93 /* The number of rows of the matrix C. M >= 0. */
94
95 /* N (input) INTEGER */
96 /* The number of columns of the matrix C. N >= 0. */
97
98 /* K (input) INTEGER */
99 /* The number of elementary reflectors whose product defines */
100 /* the matrix Q. */
101 /* If SIDE = 'L', M >= K >= 0; */
102 /* if SIDE = 'R', N >= K >= 0. */
103
104 /* L (input) INTEGER */
105 /* The number of columns of the matrix A containing */
106 /* the meaningful part of the Householder reflectors. */
107 /* If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. */
108
109 /* A (input) COMPLEX array, dimension */
110 /* (LDA,M) if SIDE = 'L', */
111 /* (LDA,N) if SIDE = 'R' */
112 /* The i-th row must contain the vector which defines the */
113 /* elementary reflector H(i), for i = 1,2,...,k, as returned by */
114 /* CTZRZF in the last k rows of its array argument A. */
115 /* A is modified by the routine but restored on exit. */
116
117 /* LDA (input) INTEGER */
118 /* The leading dimension of the array A. LDA >= max(1,K). */
119
120 /* TAU (input) COMPLEX array, dimension (K) */
121 /* TAU(i) must contain the scalar factor of the elementary */
122 /* reflector H(i), as returned by CTZRZF. */
123
124 /* C (input/output) COMPLEX array, dimension (LDC,N) */
125 /* On entry, the M-by-N matrix C. */
126 /* On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q. */
127
128 /* LDC (input) INTEGER */
129 /* The leading dimension of the array C. LDC >= max(1,M). */
130
131 /* WORK (workspace/output) COMPLEX array, dimension (LWORK) */
132 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
133
134 /* LWORK (input) INTEGER */
135 /* The dimension of the array WORK. */
136 /* If SIDE = 'L', LWORK >= max(1,N); */
137 /* if SIDE = 'R', LWORK >= max(1,M). */
138 /* For optimum performance LWORK >= N*NB if SIDE = 'L', and */
139 /* LWORK >= M*NB if SIDE = 'R', where NB is the optimal */
140 /* blocksize. */
141
142 /* If LWORK = -1, then a workspace query is assumed; the routine */
143 /* only calculates the optimal size of the WORK array, returns */
144 /* this value as the first entry of the WORK array, and no error */
145 /* message related to LWORK is issued by XERBLA. */
146
147 /* INFO (output) INTEGER */
148 /* = 0: successful exit */
149 /* < 0: if INFO = -i, the i-th argument had an illegal value */
150
151 /* Further Details */
152 /* =============== */
153
154 /* Based on contributions by */
155 /* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */
156
157 /* ===================================================================== */
158
159 /* .. Parameters .. */
160 /* .. */
161 /* .. Local Scalars .. */
162 /* .. */
163 /* .. Local Arrays .. */
164 /* .. */
165 /* .. External Functions .. */
166 /* .. */
167 /* .. External Subroutines .. */
168 /* .. */
169 /* .. Intrinsic Functions .. */
170 /* .. */
171 /* .. Executable Statements .. */
172
173 /* Test the input arguments */
174
175 /* Parameter adjustments */
176 a_dim1 = *lda;
177 a_offset = 1 + a_dim1;
178 a -= a_offset;
179 --tau;
180 c_dim1 = *ldc;
181 c_offset = 1 + c_dim1;
182 c__ -= c_offset;
183 --work;
184
185 /* Function Body */
186 *info = 0;
187 left = lsame_(side, "L", (ftnlen)1, (ftnlen)1);
188 notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
189 lquery = *lwork == -1;
190
191 /* NQ is the order of Q and NW is the minimum dimension of WORK */
192
193 if (left) {
194 nq = *m;
195 nw = *n;
196 } else {
197 nq = *n;
198 nw = *m;
199 }
200 if (! left && ! lsame_(side, "R", (ftnlen)1, (ftnlen)1)) {
201 *info = -1;
202 } else if (! notran && ! lsame_(trans, "C", (ftnlen)1, (ftnlen)1)) {
203 *info = -2;
204 } else if (*m < 0) {
205 *info = -3;
206 } else if (*n < 0) {
207 *info = -4;
208 } else if (*k < 0 || *k > nq) {
209 *info = -5;
210 } else if (*l < 0 || left && *l > *m || ! left && *l > *n) {
211 *info = -6;
212 } else if (*lda < max(1,*k)) {
213 *info = -8;
214 } else if (*ldc < max(1,*m)) {
215 *info = -11;
216 } else if (*lwork < max(1,nw) && ! lquery) {
217 *info = -13;
218 }
219
220 if (*info == 0) {
221
222 /* Determine the block size. NB may be at most NBMAX, where NBMAX */
223 /* is used to define the local array T. */
224
225 /* Computing MIN */
226 /* Writing concatenation */
227 i__3[0] = 1, a__1[0] = side;
228 i__3[1] = 1, a__1[1] = trans;
229 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
230 i__1 = 64, i__2 = ilaenv_(&c__1, "CUNMRQ", ch__1, m, n, k, &c_n1, (
231 ftnlen)6, (ftnlen)2);
232 nb = min(i__1,i__2);
233 lwkopt = max(1,nw) * nb;
234 work[1].r = (real) lwkopt, work[1].i = 0.f;
235 }
236
237 if (*info != 0) {
238 i__1 = -(*info);
239 xerbla_("CUNMRZ", &i__1, (ftnlen)6);
240 return 0;
241 } else if (lquery) {
242 return 0;
243 }
244
245 /* Quick return if possible */
246
247 if (*m == 0 || *n == 0 || *k == 0) {
248 work[1].r = 1.f, work[1].i = 0.f;
249 return 0;
250 }
251
252 /* Determine the block size. NB may be at most NBMAX, where NBMAX */
253 /* is used to define the local array T. */
254
255 /* Computing MIN */
256 /* Writing concatenation */
257 i__3[0] = 1, a__1[0] = side;
258 i__3[1] = 1, a__1[1] = trans;
259 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
260 i__1 = 64, i__2 = ilaenv_(&c__1, "CUNMRQ", ch__1, m, n, k, &c_n1, (ftnlen)
261 6, (ftnlen)2);
262 nb = min(i__1,i__2);
263 nbmin = 2;
264 ldwork = nw;
265 if (nb > 1 && nb < *k) {
266 iws = nw * nb;
267 if (*lwork < iws) {
268 nb = *lwork / ldwork;
269 /* Computing MAX */
270 /* Writing concatenation */
271 i__3[0] = 1, a__1[0] = side;
272 i__3[1] = 1, a__1[1] = trans;
273 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
274 i__1 = 2, i__2 = ilaenv_(&c__2, "CUNMRQ", ch__1, m, n, k, &c_n1, (
275 ftnlen)6, (ftnlen)2);
276 nbmin = max(i__1,i__2);
277 }
278 } else {
279 iws = nw;
280 }
281
282 if (nb < nbmin || nb >= *k) {
283
284 /* Use unblocked code */
285
286 cunmr3_(side, trans, m, n, k, l, &a[a_offset], lda, &tau[1], &c__[
287 c_offset], ldc, &work[1], &iinfo, (ftnlen)1, (ftnlen)1);
288 } else {
289
290 /* Use blocked code */
291
292 if (left && ! notran || ! left && notran) {
293 i1 = 1;
294 i2 = *k;
295 i3 = nb;
296 } else {
297 i1 = (*k - 1) / nb * nb + 1;
298 i2 = 1;
299 i3 = -nb;
300 }
301
302 if (left) {
303 ni = *n;
304 jc = 1;
305 ja = *m - *l + 1;
306 } else {
307 mi = *m;
308 ic = 1;
309 ja = *n - *l + 1;
310 }
311
312 if (notran) {
313 *(unsigned char *)transt = 'C';
314 } else {
315 *(unsigned char *)transt = 'N';
316 }
317
318 i__1 = i2;
319 i__2 = i3;
320 for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
321 /* Computing MIN */
322 i__4 = nb, i__5 = *k - i__ + 1;
323 ib = min(i__4,i__5);
324
325 /* Form the triangular factor of the block reflector */
326 /* H = H(i+ib-1) . . . H(i+1) H(i) */
327
328 clarzt_("Backward", "Rowwise", l, &ib, &a[i__ + ja * a_dim1], lda,
329 &tau[i__], t, &c__65, (ftnlen)8, (ftnlen)7);
330
331 if (left) {
332
333 /* H or H' is applied to C(i:m,1:n) */
334
335 mi = *m - i__ + 1;
336 ic = i__;
337 } else {
338
339 /* H or H' is applied to C(1:m,i:n) */
340
341 ni = *n - i__ + 1;
342 jc = i__;
343 }
344
345 /* Apply H or H' */
346
347 clarzb_(side, transt, "Backward", "Rowwise", &mi, &ni, &ib, l, &a[
348 i__ + ja * a_dim1], lda, t, &c__65, &c__[ic + jc * c_dim1]
349 , ldc, &work[1], &ldwork, (ftnlen)1, (ftnlen)1, (ftnlen)8,
350 (ftnlen)7);
351 /* L10: */
352 }
353
354 }
355
356 work[1].r = (real) lwkopt, work[1].i = 0.f;
357
358 return 0;
359
360 /* End of CUNMRZ */
361
362 } /* cunmrz_ */
363
364