1 SUBROUTINE SLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
2 *
3 * -- LAPACK auxiliary test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER INIT, SIDE
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER ISEED( 4 )
13 REAL A( LDA, * ), X( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SLAROR pre- or post-multiplies an M by N matrix A by a random
20 * orthogonal matrix U, overwriting A. A may optionally be initialized
21 * to the identity matrix before multiplying by U. U is generated using
22 * the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
23 *
24 * Arguments
25 * =========
26 *
27 * SIDE (input) CHARACTER*1
28 * Specifies whether A is multiplied on the left or right by U.
29 * = 'L': Multiply A on the left (premultiply) by U
30 * = 'R': Multiply A on the right (postmultiply) by U'
31 * = 'C' or 'T': Multiply A on the left by U and the right
32 * by U' (Here, U' means U-transpose.)
33 *
34 * INIT (input) CHARACTER*1
35 * Specifies whether or not A should be initialized to the
36 * identity matrix.
37 * = 'I': Initialize A to (a section of) the identity matrix
38 * before applying U.
39 * = 'N': No initialization. Apply U to the input matrix A.
40 *
41 * INIT = 'I' may be used to generate square or rectangular
42 * orthogonal matrices:
43 *
44 * For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
45 * to each other, as will the columns.
46 *
47 * If M < N, SIDE = 'R' produces a dense matrix whose rows are
48 * orthogonal and whose columns are not, while SIDE = 'L'
49 * produces a matrix whose rows are orthogonal, and whose first
50 * M columns are orthogonal, and whose remaining columns are
51 * zero.
52 *
53 * If M > N, SIDE = 'L' produces a dense matrix whose columns
54 * are orthogonal and whose rows are not, while SIDE = 'R'
55 * produces a matrix whose columns are orthogonal, and whose
56 * first M rows are orthogonal, and whose remaining rows are
57 * zero.
58 *
59 * M (input) INTEGER
60 * The number of rows of A.
61 *
62 * N (input) INTEGER
63 * The number of columns of A.
64 *
65 * A (input/output) REAL array, dimension (LDA, N)
66 * On entry, the array A.
67 * On exit, overwritten by U A ( if SIDE = 'L' ),
68 * or by A U ( if SIDE = 'R' ),
69 * or by U A U' ( if SIDE = 'C' or 'T').
70 *
71 * LDA (input) INTEGER
72 * The leading dimension of the array A. LDA >= max(1,M).
73 *
74 * ISEED (input/output) INTEGER array, dimension (4)
75 * On entry ISEED specifies the seed of the random number
76 * generator. The array elements should be between 0 and 4095;
77 * if not they will be reduced mod 4096. Also, ISEED(4) must
78 * be odd. The random number generator uses a linear
79 * congruential sequence limited to small integers, and so
80 * should produce machine independent random numbers. The
81 * values of ISEED are changed on exit, and can be used in the
82 * next call to SLAROR to continue the same random number
83 * sequence.
84 *
85 * X (workspace) REAL array, dimension (3*MAX( M, N ))
86 * Workspace of length
87 * 2*M + N if SIDE = 'L',
88 * 2*N + M if SIDE = 'R',
89 * 3*N if SIDE = 'C' or 'T'.
90 *
91 * INFO (output) INTEGER
92 * An error flag. It is set to:
93 * = 0: normal return
94 * < 0: if INFO = -k, the k-th argument had an illegal value
95 * = 1: if the random numbers generated by SLARND are bad.
96 *
97 * =====================================================================
98 *
99 * .. Parameters ..
100 REAL ZERO, ONE, TOOSML
101 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0,
102 $ TOOSML = 1.0E-20 )
103 * ..
104 * .. Local Scalars ..
105 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
106 REAL FACTOR, XNORM, XNORMS
107 * ..
108 * .. External Functions ..
109 LOGICAL LSAME
110 REAL SLARND, SNRM2
111 EXTERNAL LSAME, SLARND, SNRM2
112 * ..
113 * .. External Subroutines ..
114 EXTERNAL SGEMV, SGER, SLASET, SSCAL, XERBLA
115 * ..
116 * .. Intrinsic Functions ..
117 INTRINSIC ABS, SIGN
118 * ..
119 * .. Executable Statements ..
120 *
121 IF( N.EQ.0 .OR. M.EQ.0 )
122 $ RETURN
123 *
124 ITYPE = 0
125 IF( LSAME( SIDE, 'L' ) ) THEN
126 ITYPE = 1
127 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
128 ITYPE = 2
129 ELSE IF( LSAME( SIDE, 'C' ) .OR. LSAME( SIDE, 'T' ) ) THEN
130 ITYPE = 3
131 END IF
132 *
133 * Check for argument errors.
134 *
135 INFO = 0
136 IF( ITYPE.EQ.0 ) THEN
137 INFO = -1
138 ELSE IF( M.LT.0 ) THEN
139 INFO = -3
140 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
141 INFO = -4
142 ELSE IF( LDA.LT.M ) THEN
143 INFO = -6
144 END IF
145 IF( INFO.NE.0 ) THEN
146 CALL XERBLA( 'SLAROR', -INFO )
147 RETURN
148 END IF
149 *
150 IF( ITYPE.EQ.1 ) THEN
151 NXFRM = M
152 ELSE
153 NXFRM = N
154 END IF
155 *
156 * Initialize A to the identity matrix if desired
157 *
158 IF( LSAME( INIT, 'I' ) )
159 $ CALL SLASET( 'Full', M, N, ZERO, ONE, A, LDA )
160 *
161 * If no rotation possible, multiply by random +/-1
162 *
163 * Compute rotation by computing Householder transformations
164 * H(2), H(3), ..., H(nhouse)
165 *
166 DO 10 J = 1, NXFRM
167 X( J ) = ZERO
168 10 CONTINUE
169 *
170 DO 30 IXFRM = 2, NXFRM
171 KBEG = NXFRM - IXFRM + 1
172 *
173 * Generate independent normal( 0, 1 ) random numbers
174 *
175 DO 20 J = KBEG, NXFRM
176 X( J ) = SLARND( 3, ISEED )
177 20 CONTINUE
178 *
179 * Generate a Householder transformation from the random vector X
180 *
181 XNORM = SNRM2( IXFRM, X( KBEG ), 1 )
182 XNORMS = SIGN( XNORM, X( KBEG ) )
183 X( KBEG+NXFRM ) = SIGN( ONE, -X( KBEG ) )
184 FACTOR = XNORMS*( XNORMS+X( KBEG ) )
185 IF( ABS( FACTOR ).LT.TOOSML ) THEN
186 INFO = 1
187 CALL XERBLA( 'SLAROR', INFO )
188 RETURN
189 ELSE
190 FACTOR = ONE / FACTOR
191 END IF
192 X( KBEG ) = X( KBEG ) + XNORMS
193 *
194 * Apply Householder transformation to A
195 *
196 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
197 *
198 * Apply H(k) from the left.
199 *
200 CALL SGEMV( 'T', IXFRM, N, ONE, A( KBEG, 1 ), LDA,
201 $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
202 CALL SGER( IXFRM, N, -FACTOR, X( KBEG ), 1, X( 2*NXFRM+1 ),
203 $ 1, A( KBEG, 1 ), LDA )
204 *
205 END IF
206 *
207 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
208 *
209 * Apply H(k) from the right.
210 *
211 CALL SGEMV( 'N', M, IXFRM, ONE, A( 1, KBEG ), LDA,
212 $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
213 CALL SGER( M, IXFRM, -FACTOR, X( 2*NXFRM+1 ), 1, X( KBEG ),
214 $ 1, A( 1, KBEG ), LDA )
215 *
216 END IF
217 30 CONTINUE
218 *
219 X( 2*NXFRM ) = SIGN( ONE, SLARND( 3, ISEED ) )
220 *
221 * Scale the matrix A by D.
222 *
223 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
224 DO 40 IROW = 1, M
225 CALL SSCAL( N, X( NXFRM+IROW ), A( IROW, 1 ), LDA )
226 40 CONTINUE
227 END IF
228 *
229 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
230 DO 50 JCOL = 1, N
231 CALL SSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
232 50 CONTINUE
233 END IF
234 RETURN
235 *
236 * End of SLAROR
237 *
238 END
2 *
3 * -- LAPACK auxiliary test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER INIT, SIDE
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER ISEED( 4 )
13 REAL A( LDA, * ), X( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SLAROR pre- or post-multiplies an M by N matrix A by a random
20 * orthogonal matrix U, overwriting A. A may optionally be initialized
21 * to the identity matrix before multiplying by U. U is generated using
22 * the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
23 *
24 * Arguments
25 * =========
26 *
27 * SIDE (input) CHARACTER*1
28 * Specifies whether A is multiplied on the left or right by U.
29 * = 'L': Multiply A on the left (premultiply) by U
30 * = 'R': Multiply A on the right (postmultiply) by U'
31 * = 'C' or 'T': Multiply A on the left by U and the right
32 * by U' (Here, U' means U-transpose.)
33 *
34 * INIT (input) CHARACTER*1
35 * Specifies whether or not A should be initialized to the
36 * identity matrix.
37 * = 'I': Initialize A to (a section of) the identity matrix
38 * before applying U.
39 * = 'N': No initialization. Apply U to the input matrix A.
40 *
41 * INIT = 'I' may be used to generate square or rectangular
42 * orthogonal matrices:
43 *
44 * For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
45 * to each other, as will the columns.
46 *
47 * If M < N, SIDE = 'R' produces a dense matrix whose rows are
48 * orthogonal and whose columns are not, while SIDE = 'L'
49 * produces a matrix whose rows are orthogonal, and whose first
50 * M columns are orthogonal, and whose remaining columns are
51 * zero.
52 *
53 * If M > N, SIDE = 'L' produces a dense matrix whose columns
54 * are orthogonal and whose rows are not, while SIDE = 'R'
55 * produces a matrix whose columns are orthogonal, and whose
56 * first M rows are orthogonal, and whose remaining rows are
57 * zero.
58 *
59 * M (input) INTEGER
60 * The number of rows of A.
61 *
62 * N (input) INTEGER
63 * The number of columns of A.
64 *
65 * A (input/output) REAL array, dimension (LDA, N)
66 * On entry, the array A.
67 * On exit, overwritten by U A ( if SIDE = 'L' ),
68 * or by A U ( if SIDE = 'R' ),
69 * or by U A U' ( if SIDE = 'C' or 'T').
70 *
71 * LDA (input) INTEGER
72 * The leading dimension of the array A. LDA >= max(1,M).
73 *
74 * ISEED (input/output) INTEGER array, dimension (4)
75 * On entry ISEED specifies the seed of the random number
76 * generator. The array elements should be between 0 and 4095;
77 * if not they will be reduced mod 4096. Also, ISEED(4) must
78 * be odd. The random number generator uses a linear
79 * congruential sequence limited to small integers, and so
80 * should produce machine independent random numbers. The
81 * values of ISEED are changed on exit, and can be used in the
82 * next call to SLAROR to continue the same random number
83 * sequence.
84 *
85 * X (workspace) REAL array, dimension (3*MAX( M, N ))
86 * Workspace of length
87 * 2*M + N if SIDE = 'L',
88 * 2*N + M if SIDE = 'R',
89 * 3*N if SIDE = 'C' or 'T'.
90 *
91 * INFO (output) INTEGER
92 * An error flag. It is set to:
93 * = 0: normal return
94 * < 0: if INFO = -k, the k-th argument had an illegal value
95 * = 1: if the random numbers generated by SLARND are bad.
96 *
97 * =====================================================================
98 *
99 * .. Parameters ..
100 REAL ZERO, ONE, TOOSML
101 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0,
102 $ TOOSML = 1.0E-20 )
103 * ..
104 * .. Local Scalars ..
105 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
106 REAL FACTOR, XNORM, XNORMS
107 * ..
108 * .. External Functions ..
109 LOGICAL LSAME
110 REAL SLARND, SNRM2
111 EXTERNAL LSAME, SLARND, SNRM2
112 * ..
113 * .. External Subroutines ..
114 EXTERNAL SGEMV, SGER, SLASET, SSCAL, XERBLA
115 * ..
116 * .. Intrinsic Functions ..
117 INTRINSIC ABS, SIGN
118 * ..
119 * .. Executable Statements ..
120 *
121 IF( N.EQ.0 .OR. M.EQ.0 )
122 $ RETURN
123 *
124 ITYPE = 0
125 IF( LSAME( SIDE, 'L' ) ) THEN
126 ITYPE = 1
127 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
128 ITYPE = 2
129 ELSE IF( LSAME( SIDE, 'C' ) .OR. LSAME( SIDE, 'T' ) ) THEN
130 ITYPE = 3
131 END IF
132 *
133 * Check for argument errors.
134 *
135 INFO = 0
136 IF( ITYPE.EQ.0 ) THEN
137 INFO = -1
138 ELSE IF( M.LT.0 ) THEN
139 INFO = -3
140 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
141 INFO = -4
142 ELSE IF( LDA.LT.M ) THEN
143 INFO = -6
144 END IF
145 IF( INFO.NE.0 ) THEN
146 CALL XERBLA( 'SLAROR', -INFO )
147 RETURN
148 END IF
149 *
150 IF( ITYPE.EQ.1 ) THEN
151 NXFRM = M
152 ELSE
153 NXFRM = N
154 END IF
155 *
156 * Initialize A to the identity matrix if desired
157 *
158 IF( LSAME( INIT, 'I' ) )
159 $ CALL SLASET( 'Full', M, N, ZERO, ONE, A, LDA )
160 *
161 * If no rotation possible, multiply by random +/-1
162 *
163 * Compute rotation by computing Householder transformations
164 * H(2), H(3), ..., H(nhouse)
165 *
166 DO 10 J = 1, NXFRM
167 X( J ) = ZERO
168 10 CONTINUE
169 *
170 DO 30 IXFRM = 2, NXFRM
171 KBEG = NXFRM - IXFRM + 1
172 *
173 * Generate independent normal( 0, 1 ) random numbers
174 *
175 DO 20 J = KBEG, NXFRM
176 X( J ) = SLARND( 3, ISEED )
177 20 CONTINUE
178 *
179 * Generate a Householder transformation from the random vector X
180 *
181 XNORM = SNRM2( IXFRM, X( KBEG ), 1 )
182 XNORMS = SIGN( XNORM, X( KBEG ) )
183 X( KBEG+NXFRM ) = SIGN( ONE, -X( KBEG ) )
184 FACTOR = XNORMS*( XNORMS+X( KBEG ) )
185 IF( ABS( FACTOR ).LT.TOOSML ) THEN
186 INFO = 1
187 CALL XERBLA( 'SLAROR', INFO )
188 RETURN
189 ELSE
190 FACTOR = ONE / FACTOR
191 END IF
192 X( KBEG ) = X( KBEG ) + XNORMS
193 *
194 * Apply Householder transformation to A
195 *
196 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
197 *
198 * Apply H(k) from the left.
199 *
200 CALL SGEMV( 'T', IXFRM, N, ONE, A( KBEG, 1 ), LDA,
201 $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
202 CALL SGER( IXFRM, N, -FACTOR, X( KBEG ), 1, X( 2*NXFRM+1 ),
203 $ 1, A( KBEG, 1 ), LDA )
204 *
205 END IF
206 *
207 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
208 *
209 * Apply H(k) from the right.
210 *
211 CALL SGEMV( 'N', M, IXFRM, ONE, A( 1, KBEG ), LDA,
212 $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
213 CALL SGER( M, IXFRM, -FACTOR, X( 2*NXFRM+1 ), 1, X( KBEG ),
214 $ 1, A( 1, KBEG ), LDA )
215 *
216 END IF
217 30 CONTINUE
218 *
219 X( 2*NXFRM ) = SIGN( ONE, SLARND( 3, ISEED ) )
220 *
221 * Scale the matrix A by D.
222 *
223 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
224 DO 40 IROW = 1, M
225 CALL SSCAL( N, X( NXFRM+IROW ), A( IROW, 1 ), LDA )
226 40 CONTINUE
227 END IF
228 *
229 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
230 DO 50 JCOL = 1, N
231 CALL SSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
232 50 CONTINUE
233 END IF
234 RETURN
235 *
236 * End of SLAROR
237 *
238 END