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 #include <string.h>
6 #include "f2c.h"
7
8 /* Table of constant values */
9
10 static doublereal c_b9 = 0.;
11 static doublereal c_b10 = 1.;
12 static integer c__3 = 3;
13 static integer c__1 = 1;
14
dlaror_slu(char * side,char * init,integer * m,integer * n,doublereal * a,integer * lda,integer * iseed,doublereal * x,integer * info)15 /* Subroutine */ int dlaror_slu(char *side, char *init, integer *m, integer *n,
16 doublereal *a, integer *lda, integer *iseed, doublereal *x, integer *
17 info)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, i__1, i__2;
21 doublereal d__1;
22
23 /* Builtin functions */
24 double d_sign(doublereal *, doublereal *);
25
26 /* Local variables */
27 static integer kbeg;
28 extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
29 doublereal *, integer *, doublereal *, integer *, doublereal *,
30 integer *);
31 static integer jcol, irow;
32 extern doublereal dnrm2_(integer *, doublereal *, integer *);
33 static integer j;
34 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
35 integer *);
36 extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
37 doublereal *, doublereal *, integer *, doublereal *, integer *,
38 doublereal *, doublereal *, integer *);
39 static integer ixfrm, itype, nxfrm;
40 static doublereal xnorm;
41 extern doublereal dlarnd_slu(integer *, integer *);
42 extern /* Subroutine */ int dlaset_slu(char *, integer *, integer *,
43 doublereal *, doublereal *, doublereal *, integer *);
44 extern int input_error(char *, int *);
45 static doublereal factor, xnorms;
46
47
48 /* -- LAPACK auxiliary test routine (version 2.0) --
49 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
50 Courant Institute, Argonne National Lab, and Rice University
51 September 30, 1994
52
53
54 Purpose
55 =======
56
57 DLAROR pre- or post-multiplies an M by N matrix A by a random
58 orthogonal matrix U, overwriting A. A may optionally be initialized
59
60 to the identity matrix before multiplying by U. U is generated using
61
62 the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
63
64
65 Arguments
66 =========
67
68 SIDE (input) CHARACTER*1
69 Specifies whether A is multiplied on the left or right by U.
70
71 = 'L': Multiply A on the left (premultiply) by U
72 = 'R': Multiply A on the right (postmultiply) by U'
73 = 'C' or 'T': Multiply A on the left by U and the right
74 by U' (Here, U' means U-transpose.)
75
76 INIT (input) CHARACTER*1
77 Specifies whether or not A should be initialized to the
78 identity matrix.
79 = 'I': Initialize A to (a section of) the identity matrix
80 before applying U.
81 = 'N': No initialization. Apply U to the input matrix A.
82
83 INIT = 'I' may be used to generate square or rectangular
84 orthogonal matrices:
85
86 For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
87
88 to each other, as will the columns.
89
90 If M < N, SIDE = 'R' produces a dense matrix whose rows are
91 orthogonal and whose columns are not, while SIDE = 'L'
92 produces a matrix whose rows are orthogonal, and whose first
93
94 M columns are orthogonal, and whose remaining columns are
95 zero.
96
97 If M > N, SIDE = 'L' produces a dense matrix whose columns
98 are orthogonal and whose rows are not, while SIDE = 'R'
99 produces a matrix whose columns are orthogonal, and whose
100 first M rows are orthogonal, and whose remaining rows are
101 zero.
102
103 M (input) INTEGER
104 The number of rows of A.
105
106 N (input) INTEGER
107 The number of columns of A.
108
109 A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
110 On entry, the array A.
111 On exit, overwritten by U A ( if SIDE = 'L' ),
112 or by A U ( if SIDE = 'R' ),
113 or by U A U' ( if SIDE = 'C' or 'T').
114
115 LDA (input) INTEGER
116 The leading dimension of the array A. LDA >= max(1,M).
117
118 ISEED (input/output) INTEGER array, dimension (4)
119 On entry ISEED specifies the seed of the random number
120 generator. The array elements should be between 0 and 4095;
121 if not they will be reduced mod 4096. Also, ISEED(4) must
122 be odd. The random number generator uses a linear
123 congruential sequence limited to small integers, and so
124 should produce machine independent random numbers. The
125 values of ISEED are changed on exit, and can be used in the
126 next call to DLAROR to continue the same random number
127 sequence.
128
129 X (workspace) DOUBLE PRECISION array, dimension (3*MAX( M, N ))
130
131 Workspace of length
132 2*M + N if SIDE = 'L',
133 2*N + M if SIDE = 'R',
134 3*N if SIDE = 'C' or 'T'.
135
136 INFO (output) INTEGER
137 An error flag. It is set to:
138 = 0: normal return
139 < 0: if INFO = -k, the k-th argument had an illegal value
140 = 1: if the random numbers generated by DLARND are bad.
141
142 =====================================================================
143
144
145
146 Parameter adjustments */
147 a_dim1 = *lda;
148 a_offset = a_dim1 + 1;
149 a -= a_offset;
150 --iseed;
151 --x;
152
153 /* Function Body */
154 if (*n == 0 || *m == 0) {
155 return 0;
156 }
157
158 itype = 0;
159 if (strncmp(side, "L", 1)==0) {
160 itype = 1;
161 } else if (strncmp(side, "R", 1)==0) {
162 itype = 2;
163 } else if (strncmp(side, "C", 1)==0 || strncmp(side, "T", 1)==0) {
164 itype = 3;
165 }
166
167 /* Check for argument errors. */
168
169 *info = 0;
170 if (itype == 0) {
171 *info = -1;
172 } else if (*m < 0) {
173 *info = -3;
174 } else if (*n < 0 || itype == 3 && *n != *m) {
175 *info = -4;
176 } else if (*lda < *m) {
177 *info = -6;
178 }
179 if (*info != 0) {
180 i__1 = -(*info);
181 input_error("DLAROR", &i__1);
182 return 0;
183 }
184
185 if (itype == 1) {
186 nxfrm = *m;
187 } else {
188 nxfrm = *n;
189 }
190
191 /* Initialize A to the identity matrix if desired */
192
193 if (strncmp(init, "I", 1)==0) {
194 dlaset_slu("Full", m, n, &c_b9, &c_b10, &a[a_offset], lda);
195 }
196
197 /* If no rotation possible, multiply by random +/-1
198
199 Compute rotation by computing Householder transformations
200 H(2), H(3), ..., H(nhouse) */
201
202 i__1 = nxfrm;
203 for (j = 1; j <= i__1; ++j) {
204 x[j] = 0.;
205 /* L10: */
206 }
207
208 i__1 = nxfrm;
209 for (ixfrm = 2; ixfrm <= i__1; ++ixfrm) {
210 kbeg = nxfrm - ixfrm + 1;
211
212 /* Generate independent normal( 0, 1 ) random numbers */
213
214 i__2 = nxfrm;
215 for (j = kbeg; j <= i__2; ++j) {
216 x[j] = dlarnd_slu(&c__3, &iseed[1]);
217 /* L20: */
218 }
219
220 /* Generate a Householder transformation from the random vector
221 X */
222
223 xnorm = dnrm2_(&ixfrm, &x[kbeg], &c__1);
224 xnorms = d_sign(&xnorm, &x[kbeg]);
225 d__1 = -x[kbeg];
226 x[kbeg + nxfrm] = d_sign(&c_b10, &d__1);
227 factor = xnorms * (xnorms + x[kbeg]);
228 if (abs(factor) < 1e-20) {
229 *info = 1;
230 input_error("DLAROR", info);
231 return 0;
232 } else {
233 factor = 1. / factor;
234 }
235 x[kbeg] += xnorms;
236
237 /* Apply Householder transformation to A */
238
239 if (itype == 1 || itype == 3) {
240
241 /* Apply H(k) from the left. */
242
243 dgemv_("T", &ixfrm, n, &c_b10, &a[kbeg + a_dim1], lda, &x[kbeg], &
244 c__1, &c_b9, &x[(nxfrm << 1) + 1], &c__1);
245 d__1 = -factor;
246 dger_(&ixfrm, n, &d__1, &x[kbeg], &c__1, &x[(nxfrm << 1) + 1], &
247 c__1, &a[kbeg + a_dim1], lda);
248
249 }
250
251 if (itype == 2 || itype == 3) {
252
253 /* Apply H(k) from the right. */
254
255 dgemv_("N", m, &ixfrm, &c_b10, &a[kbeg * a_dim1 + 1], lda, &x[
256 kbeg], &c__1, &c_b9, &x[(nxfrm << 1) + 1], &c__1);
257 d__1 = -factor;
258 dger_(m, &ixfrm, &d__1, &x[(nxfrm << 1) + 1], &c__1, &x[kbeg], &
259 c__1, &a[kbeg * a_dim1 + 1], lda);
260
261 }
262 /* L30: */
263 }
264
265 d__1 = dlarnd_slu(&c__3, &iseed[1]);
266 x[nxfrm * 2] = d_sign(&c_b10, &d__1);
267
268 /* Scale the matrix A by D. */
269
270 if (itype == 1 || itype == 3) {
271 i__1 = *m;
272 for (irow = 1; irow <= i__1; ++irow) {
273 dscal_(n, &x[nxfrm + irow], &a[irow + a_dim1], lda);
274 /* L40: */
275 }
276 }
277
278 if (itype == 2 || itype == 3) {
279 i__1 = *n;
280 for (jcol = 1; jcol <= i__1; ++jcol) {
281 dscal_(m, &x[nxfrm + jcol], &a[jcol * a_dim1 + 1], &c__1);
282 /* L50: */
283 }
284 }
285 return 0;
286
287 /* End of DLAROR */
288
289 } /* dlaror_slu */
290
291