1 /* -- translated by f2c (version 19940927).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.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__3 = 3;
13 static integer c__1 = 1;
14
claror_(char * side,char * init,integer * m,integer * n,complex * a,integer * lda,integer * iseed,complex * x,integer * info)15 /* Subroutine */ int claror_(char *side, char *init, integer *m, integer *n,
16 complex *a, integer *lda, integer *iseed, complex *x, integer *info)
17 {
18 /* System generated locals */
19 integer a_dim1, a_offset, i__1, i__2, i__3;
20 complex q__1, q__2;
21
22 /* Builtin functions */
23 double c_abs(complex *);
24 void r_cnjg(complex *, complex *);
25
26 /* Local variables */
27 static integer kbeg, jcol;
28 static real xabs;
29 static integer irow, j;
30 extern /* Subroutine */ int cgerc_(integer *, integer *, complex *,
31 complex *, integer *, complex *, integer *, complex *, integer *),
32 cscal_(integer *, complex *, complex *, integer *);
33 extern logical lsame_(char *, char *);
34 extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
35 , complex *, integer *, complex *, integer *, complex *, complex *
36 , integer *);
37 static complex csign;
38 static integer ixfrm, itype, nxfrm;
39 static real xnorm;
40 extern real scnrm2_(integer *, complex *, integer *);
41 extern /* Subroutine */ int clacgv_(integer *, complex *, integer *);
42 extern /* Complex */ VOID clarnd_(complex *, integer *, integer *);
43 extern /* Subroutine */ int claset_(char *, integer *, integer *, complex
44 *, complex *, complex *, integer *), xerbla_(char *,
45 integer *);
46 static real factor;
47 static complex xnorms;
48
49
50 /* -- LAPACK auxiliary test routine (version 2.0) --
51 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
52 Courant Institute, Argonne National Lab, and Rice University
53 September 30, 1994
54
55
56 Purpose
57 =======
58
59 CLAROR pre- or post-multiplies an M by N matrix A by a random
60 unitary matrix U, overwriting A. A may optionally be
61 initialized to the identity matrix before multiplying by U.
62 U is generated using the method of G.W. Stewart
63 ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ).
64 (BLAS-2 version)
65
66 Arguments
67 =========
68
69 SIDE - CHARACTER*1
70 SIDE specifies whether A is multiplied on the left or right
71
72 by U.
73 SIDE = 'L' Multiply A on the left (premultiply) by U
74 SIDE = 'R' Multiply A on the right (postmultiply) by U*
75 SIDE = 'C' Multiply A on the left by U and the right by U*
76 SIDE = 'T' Multiply A on the left by U and the right by U'
77 Not modified.
78
79 INIT - CHARACTER*1
80 INIT specifies whether or not A should be initialized to
81 the identity matrix.
82 INIT = 'I' Initialize A to (a section of) the
83 identity matrix before applying U.
84 INIT = 'N' No initialization. Apply U to the
85 input matrix A.
86
87 INIT = 'I' may be used to generate square (i.e., unitary)
88 or rectangular orthogonal matrices (orthogonality being
89 in the sense of CDOTC):
90
91 For square matrices, M=N, and SIDE many be either 'L' or
92 'R'; the rows will be orthogonal to each other, as will the
93
94 columns.
95 For rectangular matrices where M < N, SIDE = 'R' will
96 produce a dense matrix whose rows will be orthogonal and
97 whose columns will not, while SIDE = 'L' will produce a
98 matrix whose rows will be orthogonal, and whose first M
99 columns will be orthogonal, the remaining columns being
100 zero.
101 For matrices where M > N, just use the previous
102 explaination, interchanging 'L' and 'R' and "rows" and
103 "columns".
104
105 Not modified.
106
107 M - INTEGER
108 Number of rows of A. Not modified.
109
110 N - INTEGER
111 Number of columns of A. Not modified.
112
113 A - COMPLEX array, dimension ( LDA, N )
114 Input and output array. Overwritten by U A ( if SIDE = 'L' )
115
116 or by A U ( if SIDE = 'R' )
117 or by U A U* ( if SIDE = 'C')
118 or by U A U' ( if SIDE = 'T') on exit.
119
120 LDA - INTEGER
121 Leading dimension of A. Must be at least MAX ( 1, M ).
122 Not modified.
123
124 ISEED - INTEGER array, dimension ( 4 )
125 On entry ISEED specifies the seed of the random number
126 generator. The array elements should be between 0 and 4095;
127
128 if not they will be reduced mod 4096. Also, ISEED(4) must
129 be odd. The random number generator uses a linear
130 congruential sequence limited to small integers, and so
131 should produce machine independent random numbers. The
132 values of ISEED are changed on exit, and can be used in the
133
134 next call to CLAROR to continue the same random number
135 sequence.
136 Modified.
137
138 X - COMPLEX array, dimension ( 3*MAX( M, N ) )
139 Workspace. Of length:
140 2*M + N if SIDE = 'L',
141 2*N + M if SIDE = 'R',
142 3*N if SIDE = 'C' or 'T'.
143 Modified.
144
145 INFO - INTEGER
146 An error flag. It is set to:
147 0 if no error.
148 1 if CLARND returned a bad random number (installation
149 problem)
150 -1 if SIDE is not L, R, C, or T.
151 -3 if M is negative.
152 -4 if N is negative or if SIDE is C or T and N is not equal
153
154 to M.
155 -6 if LDA is less than M.
156
157 =====================================================================
158
159
160
161 Parameter adjustments */
162 a_dim1 = *lda;
163 a_offset = a_dim1 + 1;
164 a -= a_offset;
165 --iseed;
166 --x;
167
168 /* Function Body */
169 if (*n == 0 || *m == 0) {
170 return 0;
171 }
172
173 itype = 0;
174 if (lsame_(side, "L")) {
175 itype = 1;
176 } else if (lsame_(side, "R")) {
177 itype = 2;
178 } else if (lsame_(side, "C")) {
179 itype = 3;
180 } else if (lsame_(side, "T")) {
181 itype = 4;
182 }
183
184 /* Check for argument errors. */
185
186 *info = 0;
187 if (itype == 0) {
188 *info = -1;
189 } else if (*m < 0) {
190 *info = -3;
191 } else if (*n < 0 || itype == 3 && *n != *m) {
192 *info = -4;
193 } else if (*lda < *m) {
194 *info = -6;
195 }
196 if (*info != 0) {
197 i__1 = -(*info);
198 xerbla_("CLAROR", &i__1);
199 return 0;
200 }
201
202 if (itype == 1) {
203 nxfrm = *m;
204 } else {
205 nxfrm = *n;
206 }
207
208 /* Initialize A to the identity matrix if desired */
209
210 if (lsame_(init, "I")) {
211 claset_("Full", m, n, &c_b1, &c_b2, &a[a_offset], lda);
212 }
213
214 /* If no rotation possible, still multiply by
215 a random complex number from the circle |x| = 1
216
217 2) Compute Rotation by computing Householder
218 Transformations H(2), H(3), ..., H(n). Note that the
219 order in which they are computed is irrelevant. */
220
221 i__1 = nxfrm;
222 for (j = 1; j <= i__1; ++j) {
223 i__2 = j;
224 x[i__2].r = 0.f, x[i__2].i = 0.f;
225 /* L40: */
226 }
227
228 i__1 = nxfrm;
229 for (ixfrm = 2; ixfrm <= i__1; ++ixfrm) {
230 kbeg = nxfrm - ixfrm + 1;
231
232 /* Generate independent normal( 0, 1 ) random numbers */
233
234 i__2 = nxfrm;
235 for (j = kbeg; j <= i__2; ++j) {
236 i__3 = j;
237 clarnd_(&q__1, &c__3, &iseed[1]);
238 x[i__3].r = q__1.r, x[i__3].i = q__1.i;
239 /* L50: */
240 }
241
242 /* Generate a Householder transformation from the random vector
243 X */
244
245 xnorm = scnrm2_(&ixfrm, &x[kbeg], &c__1);
246 xabs = c_abs(&x[kbeg]);
247 if (xabs != 0.f) {
248 i__2 = kbeg;
249 q__1.r = x[i__2].r / xabs, q__1.i = x[i__2].i / xabs;
250 csign.r = q__1.r, csign.i = q__1.i;
251 } else {
252 csign.r = 1.f, csign.i = 0.f;
253 }
254 q__1.r = xnorm * csign.r, q__1.i = xnorm * csign.i;
255 xnorms.r = q__1.r, xnorms.i = q__1.i;
256 i__2 = nxfrm + kbeg;
257 q__1.r = -(doublereal)csign.r, q__1.i = -(doublereal)csign.i;
258 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
259 factor = xnorm * (xnorm + xabs);
260 if (dabs(factor) < 1e-20f) {
261 *info = 1;
262 i__2 = -(*info);
263 xerbla_("CLAROR", &i__2);
264 return 0;
265 } else {
266 factor = 1.f / factor;
267 }
268 i__2 = kbeg;
269 i__3 = kbeg;
270 q__1.r = x[i__3].r + xnorms.r, q__1.i = x[i__3].i + xnorms.i;
271 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
272
273 /* Apply Householder transformation to A */
274
275 if (itype == 1 || itype == 3 || itype == 4) {
276
277 /* Apply H(k) on the left of A */
278
279 cgemv_("C", &ixfrm, n, &c_b2, &a[kbeg + a_dim1], lda, &x[kbeg], &
280 c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
281 q__2.r = factor, q__2.i = 0.f;
282 q__1.r = -(doublereal)q__2.r, q__1.i = -(doublereal)q__2.i;
283 cgerc_(&ixfrm, n, &q__1, &x[kbeg], &c__1, &x[(nxfrm << 1) + 1], &
284 c__1, &a[kbeg + a_dim1], lda);
285
286 }
287
288 if (itype >= 2 && itype <= 4) {
289
290 /* Apply H(k)* (or H(k)') on the right of A */
291
292 if (itype == 4) {
293 clacgv_(&ixfrm, &x[kbeg], &c__1);
294 }
295
296 cgemv_("N", m, &ixfrm, &c_b2, &a[kbeg * a_dim1 + 1], lda, &x[kbeg]
297 , &c__1, &c_b1, &x[(nxfrm << 1) + 1], &c__1);
298 q__2.r = factor, q__2.i = 0.f;
299 q__1.r = -(doublereal)q__2.r, q__1.i = -(doublereal)q__2.i;
300 cgerc_(m, &ixfrm, &q__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
301 c__1, &a[kbeg * a_dim1 + 1], lda);
302
303 }
304 /* L60: */
305 }
306
307 clarnd_(&q__1, &c__3, &iseed[1]);
308 x[1].r = q__1.r, x[1].i = q__1.i;
309 xabs = c_abs(&x[1]);
310 if (xabs != 0.f) {
311 q__1.r = x[1].r / xabs, q__1.i = x[1].i / xabs;
312 csign.r = q__1.r, csign.i = q__1.i;
313 } else {
314 csign.r = 1.f, csign.i = 0.f;
315 }
316 i__1 = nxfrm << 1;
317 x[i__1].r = csign.r, x[i__1].i = csign.i;
318
319 /* Scale the matrix A by D. */
320
321 if (itype == 1 || itype == 3 || itype == 4) {
322 i__1 = *m;
323 for (irow = 1; irow <= i__1; ++irow) {
324 r_cnjg(&q__1, &x[nxfrm + irow]);
325 cscal_(n, &q__1, &a[irow + a_dim1], lda);
326 /* L70: */
327 }
328 }
329
330 if (itype == 2 || itype == 3) {
331 i__1 = *n;
332 for (jcol = 1; jcol <= i__1; ++jcol) {
333 cscal_(m, &x[nxfrm + jcol], &a[jcol * a_dim1 + 1], &c__1);
334 /* L80: */
335 }
336 }
337
338 if (itype == 4) {
339 i__1 = *n;
340 for (jcol = 1; jcol <= i__1; ++jcol) {
341 r_cnjg(&q__1, &x[nxfrm + jcol]);
342 cscal_(m, &q__1, &a[jcol * a_dim1 + 1], &c__1);
343 /* L90: */
344 }
345 }
346 return 0;
347
348 /* End of CLAROR */
349
350 } /* claror_ */
351
352