1 SUBROUTINE ZLAROR( 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 * June 2010
6 *
7 * .. Scalar Arguments ..
8 CHARACTER INIT, SIDE
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER ISEED( 4 )
13 COMPLEX*16 A( LDA, * ), X( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLAROR pre- or post-multiplies an M by N matrix A by a random
20 * unitary matrix U, overwriting A. A may optionally be
21 * initialized to the identity matrix before multiplying by U.
22 * U is generated using the method of G.W. Stewart
23 * ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ).
24 * (BLAS-2 version)
25 *
26 * Arguments
27 * =========
28 *
29 * SIDE (input) CHARACTER*1
30 * SIDE specifies whether A is multiplied on the left or right
31 * by U.
32 * SIDE = 'L' Multiply A on the left (premultiply) by U
33 * SIDE = 'R' Multiply A on the right (postmultiply) by U*
34 * SIDE = 'C' Multiply A on the left by U and the right by U*
35 * SIDE = 'T' Multiply A on the left by U and the right by U'
36 * Not modified.
37 *
38 * INIT (input) CHARACTER*1
39 * INIT specifies whether or not A should be initialized to
40 * the identity matrix.
41 * INIT = 'I' Initialize A to (a section of) the
42 * identity matrix before applying U.
43 * INIT = 'N' No initialization. Apply U to the
44 * input matrix A.
45 *
46 * INIT = 'I' may be used to generate square (i.e., unitary)
47 * or rectangular orthogonal matrices (orthogonality being
48 * in the sense of ZDOTC):
49 *
50 * For square matrices, M=N, and SIDE many be either 'L' or
51 * 'R'; the rows will be orthogonal to each other, as will the
52 * columns.
53 * For rectangular matrices where M < N, SIDE = 'R' will
54 * produce a dense matrix whose rows will be orthogonal and
55 * whose columns will not, while SIDE = 'L' will produce a
56 * matrix whose rows will be orthogonal, and whose first M
57 * columns will be orthogonal, the remaining columns being
58 * zero.
59 * For matrices where M > N, just use the previous
60 * explanation, interchanging 'L' and 'R' and "rows" and
61 * "columns".
62 *
63 * Not modified.
64 *
65 * M (input) INTEGER
66 * Number of rows of A. Not modified.
67 *
68 * N (input) INTEGER
69 * Number of columns of A. Not modified.
70 *
71 * A COMPLEX*16 array, dimension ( LDA, N )
72 * Input and output array. Overwritten by U A ( if SIDE = 'L' )
73 * or by A U ( if SIDE = 'R' )
74 * or by U A U* ( if SIDE = 'C')
75 * or by U A U' ( if SIDE = 'T') on exit.
76 *
77 * LDA (input) INTEGER
78 * Leading dimension of A. Must be at least MAX ( 1, M ).
79 * Not modified.
80 *
81 * ISEED (input/output) INTEGER array, dimension ( 4 )
82 * On entry ISEED specifies the seed of the random number
83 * generator. The array elements should be between 0 and 4095;
84 * if not they will be reduced mod 4096. Also, ISEED(4) must
85 * be odd. The random number generator uses a linear
86 * congruential sequence limited to small integers, and so
87 * should produce machine independent random numbers. The
88 * values of ISEED are changed on exit, and can be used in the
89 * next call to ZLAROR to continue the same random number
90 * sequence.
91 * Modified.
92 *
93 * X (workspace) COMPLEX*16 array, dimension ( 3*MAX( M, N ) )
94 * Workspace. Of length:
95 * 2*M + N if SIDE = 'L',
96 * 2*N + M if SIDE = 'R',
97 * 3*N if SIDE = 'C' or 'T'.
98 * Modified.
99 *
100 * INFO (output) INTEGER
101 * An error flag. It is set to:
102 * 0 if no error.
103 * 1 if ZLARND returned a bad random number (installation
104 * problem)
105 * -1 if SIDE is not L, R, C, or T.
106 * -3 if M is negative.
107 * -4 if N is negative or if SIDE is C or T and N is not equal
108 * to M.
109 * -6 if LDA is less than M.
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114 DOUBLE PRECISION ZERO, ONE, TOOSML
115 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
116 $ TOOSML = 1.0D-20 )
117 COMPLEX*16 CZERO, CONE
118 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
119 $ CONE = ( 1.0D+0, 0.0D+0 ) )
120 * ..
121 * .. Local Scalars ..
122 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
123 DOUBLE PRECISION FACTOR, XABS, XNORM
124 COMPLEX*16 CSIGN, XNORMS
125 * ..
126 * .. External Functions ..
127 LOGICAL LSAME
128 DOUBLE PRECISION DZNRM2
129 COMPLEX*16 ZLARND
130 EXTERNAL LSAME, DZNRM2, ZLARND
131 * ..
132 * .. External Subroutines ..
133 EXTERNAL XERBLA, ZGEMV, ZGERC, ZLACGV, ZLASET, ZSCAL
134 * ..
135 * .. Intrinsic Functions ..
136 INTRINSIC ABS, DCMPLX, DCONJG
137 * ..
138 * .. Executable Statements ..
139 *
140 IF( N.EQ.0 .OR. M.EQ.0 )
141 $ RETURN
142 *
143 ITYPE = 0
144 IF( LSAME( SIDE, 'L' ) ) THEN
145 ITYPE = 1
146 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
147 ITYPE = 2
148 ELSE IF( LSAME( SIDE, 'C' ) ) THEN
149 ITYPE = 3
150 ELSE IF( LSAME( SIDE, 'T' ) ) THEN
151 ITYPE = 4
152 END IF
153 *
154 * Check for argument errors.
155 *
156 INFO = 0
157 IF( ITYPE.EQ.0 ) THEN
158 INFO = -1
159 ELSE IF( M.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
162 INFO = -4
163 ELSE IF( LDA.LT.M ) THEN
164 INFO = -6
165 END IF
166 IF( INFO.NE.0 ) THEN
167 CALL XERBLA( 'ZLAROR', -INFO )
168 RETURN
169 END IF
170 *
171 IF( ITYPE.EQ.1 ) THEN
172 NXFRM = M
173 ELSE
174 NXFRM = N
175 END IF
176 *
177 * Initialize A to the identity matrix if desired
178 *
179 IF( LSAME( INIT, 'I' ) )
180 $ CALL ZLASET( 'Full', M, N, CZERO, CONE, A, LDA )
181 *
182 * If no rotation possible, still multiply by
183 * a random complex number from the circle |x| = 1
184 *
185 * 2) Compute Rotation by computing Householder
186 * Transformations H(2), H(3), ..., H(n). Note that the
187 * order in which they are computed is irrelevant.
188 *
189 DO 10 J = 1, NXFRM
190 X( J ) = CZERO
191 10 CONTINUE
192 *
193 DO 30 IXFRM = 2, NXFRM
194 KBEG = NXFRM - IXFRM + 1
195 *
196 * Generate independent normal( 0, 1 ) random numbers
197 *
198 DO 20 J = KBEG, NXFRM
199 X( J ) = ZLARND( 3, ISEED )
200 20 CONTINUE
201 *
202 * Generate a Householder transformation from the random vector X
203 *
204 XNORM = DZNRM2( IXFRM, X( KBEG ), 1 )
205 XABS = ABS( X( KBEG ) )
206 IF( XABS.NE.CZERO ) THEN
207 CSIGN = X( KBEG ) / XABS
208 ELSE
209 CSIGN = CONE
210 END IF
211 XNORMS = CSIGN*XNORM
212 X( NXFRM+KBEG ) = -CSIGN
213 FACTOR = XNORM*( XNORM+XABS )
214 IF( ABS( FACTOR ).LT.TOOSML ) THEN
215 INFO = 1
216 CALL XERBLA( 'ZLAROR', -INFO )
217 RETURN
218 ELSE
219 FACTOR = ONE / FACTOR
220 END IF
221 X( KBEG ) = X( KBEG ) + XNORMS
222 *
223 * Apply Householder transformation to A
224 *
225 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN
226 *
227 * Apply H(k) on the left of A
228 *
229 CALL ZGEMV( 'C', IXFRM, N, CONE, A( KBEG, 1 ), LDA,
230 $ X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 )
231 CALL ZGERC( IXFRM, N, -DCMPLX( FACTOR ), X( KBEG ), 1,
232 $ X( 2*NXFRM+1 ), 1, A( KBEG, 1 ), LDA )
233 *
234 END IF
235 *
236 IF( ITYPE.GE.2 .AND. ITYPE.LE.4 ) THEN
237 *
238 * Apply H(k)* (or H(k)') on the right of A
239 *
240 IF( ITYPE.EQ.4 ) THEN
241 CALL ZLACGV( IXFRM, X( KBEG ), 1 )
242 END IF
243 *
244 CALL ZGEMV( 'N', M, IXFRM, CONE, A( 1, KBEG ), LDA,
245 $ X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 )
246 CALL ZGERC( M, IXFRM, -DCMPLX( FACTOR ), X( 2*NXFRM+1 ), 1,
247 $ X( KBEG ), 1, A( 1, KBEG ), LDA )
248 *
249 END IF
250 30 CONTINUE
251 *
252 X( 1 ) = ZLARND( 3, ISEED )
253 XABS = ABS( X( 1 ) )
254 IF( XABS.NE.ZERO ) THEN
255 CSIGN = X( 1 ) / XABS
256 ELSE
257 CSIGN = CONE
258 END IF
259 X( 2*NXFRM ) = CSIGN
260 *
261 * Scale the matrix A by D.
262 *
263 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN
264 DO 40 IROW = 1, M
265 CALL ZSCAL( N, DCONJG( X( NXFRM+IROW ) ), A( IROW, 1 ),
266 $ LDA )
267 40 CONTINUE
268 END IF
269 *
270 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
271 DO 50 JCOL = 1, N
272 CALL ZSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
273 50 CONTINUE
274 END IF
275 *
276 IF( ITYPE.EQ.4 ) THEN
277 DO 60 JCOL = 1, N
278 CALL ZSCAL( M, DCONJG( X( NXFRM+JCOL ) ), A( 1, JCOL ), 1 )
279 60 CONTINUE
280 END IF
281 RETURN
282 *
283 * End of ZLAROR
284 *
285 END
2 *
3 * -- LAPACK auxiliary test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * June 2010
6 *
7 * .. Scalar Arguments ..
8 CHARACTER INIT, SIDE
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER ISEED( 4 )
13 COMPLEX*16 A( LDA, * ), X( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLAROR pre- or post-multiplies an M by N matrix A by a random
20 * unitary matrix U, overwriting A. A may optionally be
21 * initialized to the identity matrix before multiplying by U.
22 * U is generated using the method of G.W. Stewart
23 * ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ).
24 * (BLAS-2 version)
25 *
26 * Arguments
27 * =========
28 *
29 * SIDE (input) CHARACTER*1
30 * SIDE specifies whether A is multiplied on the left or right
31 * by U.
32 * SIDE = 'L' Multiply A on the left (premultiply) by U
33 * SIDE = 'R' Multiply A on the right (postmultiply) by U*
34 * SIDE = 'C' Multiply A on the left by U and the right by U*
35 * SIDE = 'T' Multiply A on the left by U and the right by U'
36 * Not modified.
37 *
38 * INIT (input) CHARACTER*1
39 * INIT specifies whether or not A should be initialized to
40 * the identity matrix.
41 * INIT = 'I' Initialize A to (a section of) the
42 * identity matrix before applying U.
43 * INIT = 'N' No initialization. Apply U to the
44 * input matrix A.
45 *
46 * INIT = 'I' may be used to generate square (i.e., unitary)
47 * or rectangular orthogonal matrices (orthogonality being
48 * in the sense of ZDOTC):
49 *
50 * For square matrices, M=N, and SIDE many be either 'L' or
51 * 'R'; the rows will be orthogonal to each other, as will the
52 * columns.
53 * For rectangular matrices where M < N, SIDE = 'R' will
54 * produce a dense matrix whose rows will be orthogonal and
55 * whose columns will not, while SIDE = 'L' will produce a
56 * matrix whose rows will be orthogonal, and whose first M
57 * columns will be orthogonal, the remaining columns being
58 * zero.
59 * For matrices where M > N, just use the previous
60 * explanation, interchanging 'L' and 'R' and "rows" and
61 * "columns".
62 *
63 * Not modified.
64 *
65 * M (input) INTEGER
66 * Number of rows of A. Not modified.
67 *
68 * N (input) INTEGER
69 * Number of columns of A. Not modified.
70 *
71 * A COMPLEX*16 array, dimension ( LDA, N )
72 * Input and output array. Overwritten by U A ( if SIDE = 'L' )
73 * or by A U ( if SIDE = 'R' )
74 * or by U A U* ( if SIDE = 'C')
75 * or by U A U' ( if SIDE = 'T') on exit.
76 *
77 * LDA (input) INTEGER
78 * Leading dimension of A. Must be at least MAX ( 1, M ).
79 * Not modified.
80 *
81 * ISEED (input/output) INTEGER array, dimension ( 4 )
82 * On entry ISEED specifies the seed of the random number
83 * generator. The array elements should be between 0 and 4095;
84 * if not they will be reduced mod 4096. Also, ISEED(4) must
85 * be odd. The random number generator uses a linear
86 * congruential sequence limited to small integers, and so
87 * should produce machine independent random numbers. The
88 * values of ISEED are changed on exit, and can be used in the
89 * next call to ZLAROR to continue the same random number
90 * sequence.
91 * Modified.
92 *
93 * X (workspace) COMPLEX*16 array, dimension ( 3*MAX( M, N ) )
94 * Workspace. Of length:
95 * 2*M + N if SIDE = 'L',
96 * 2*N + M if SIDE = 'R',
97 * 3*N if SIDE = 'C' or 'T'.
98 * Modified.
99 *
100 * INFO (output) INTEGER
101 * An error flag. It is set to:
102 * 0 if no error.
103 * 1 if ZLARND returned a bad random number (installation
104 * problem)
105 * -1 if SIDE is not L, R, C, or T.
106 * -3 if M is negative.
107 * -4 if N is negative or if SIDE is C or T and N is not equal
108 * to M.
109 * -6 if LDA is less than M.
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114 DOUBLE PRECISION ZERO, ONE, TOOSML
115 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
116 $ TOOSML = 1.0D-20 )
117 COMPLEX*16 CZERO, CONE
118 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
119 $ CONE = ( 1.0D+0, 0.0D+0 ) )
120 * ..
121 * .. Local Scalars ..
122 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
123 DOUBLE PRECISION FACTOR, XABS, XNORM
124 COMPLEX*16 CSIGN, XNORMS
125 * ..
126 * .. External Functions ..
127 LOGICAL LSAME
128 DOUBLE PRECISION DZNRM2
129 COMPLEX*16 ZLARND
130 EXTERNAL LSAME, DZNRM2, ZLARND
131 * ..
132 * .. External Subroutines ..
133 EXTERNAL XERBLA, ZGEMV, ZGERC, ZLACGV, ZLASET, ZSCAL
134 * ..
135 * .. Intrinsic Functions ..
136 INTRINSIC ABS, DCMPLX, DCONJG
137 * ..
138 * .. Executable Statements ..
139 *
140 IF( N.EQ.0 .OR. M.EQ.0 )
141 $ RETURN
142 *
143 ITYPE = 0
144 IF( LSAME( SIDE, 'L' ) ) THEN
145 ITYPE = 1
146 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
147 ITYPE = 2
148 ELSE IF( LSAME( SIDE, 'C' ) ) THEN
149 ITYPE = 3
150 ELSE IF( LSAME( SIDE, 'T' ) ) THEN
151 ITYPE = 4
152 END IF
153 *
154 * Check for argument errors.
155 *
156 INFO = 0
157 IF( ITYPE.EQ.0 ) THEN
158 INFO = -1
159 ELSE IF( M.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
162 INFO = -4
163 ELSE IF( LDA.LT.M ) THEN
164 INFO = -6
165 END IF
166 IF( INFO.NE.0 ) THEN
167 CALL XERBLA( 'ZLAROR', -INFO )
168 RETURN
169 END IF
170 *
171 IF( ITYPE.EQ.1 ) THEN
172 NXFRM = M
173 ELSE
174 NXFRM = N
175 END IF
176 *
177 * Initialize A to the identity matrix if desired
178 *
179 IF( LSAME( INIT, 'I' ) )
180 $ CALL ZLASET( 'Full', M, N, CZERO, CONE, A, LDA )
181 *
182 * If no rotation possible, still multiply by
183 * a random complex number from the circle |x| = 1
184 *
185 * 2) Compute Rotation by computing Householder
186 * Transformations H(2), H(3), ..., H(n). Note that the
187 * order in which they are computed is irrelevant.
188 *
189 DO 10 J = 1, NXFRM
190 X( J ) = CZERO
191 10 CONTINUE
192 *
193 DO 30 IXFRM = 2, NXFRM
194 KBEG = NXFRM - IXFRM + 1
195 *
196 * Generate independent normal( 0, 1 ) random numbers
197 *
198 DO 20 J = KBEG, NXFRM
199 X( J ) = ZLARND( 3, ISEED )
200 20 CONTINUE
201 *
202 * Generate a Householder transformation from the random vector X
203 *
204 XNORM = DZNRM2( IXFRM, X( KBEG ), 1 )
205 XABS = ABS( X( KBEG ) )
206 IF( XABS.NE.CZERO ) THEN
207 CSIGN = X( KBEG ) / XABS
208 ELSE
209 CSIGN = CONE
210 END IF
211 XNORMS = CSIGN*XNORM
212 X( NXFRM+KBEG ) = -CSIGN
213 FACTOR = XNORM*( XNORM+XABS )
214 IF( ABS( FACTOR ).LT.TOOSML ) THEN
215 INFO = 1
216 CALL XERBLA( 'ZLAROR', -INFO )
217 RETURN
218 ELSE
219 FACTOR = ONE / FACTOR
220 END IF
221 X( KBEG ) = X( KBEG ) + XNORMS
222 *
223 * Apply Householder transformation to A
224 *
225 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN
226 *
227 * Apply H(k) on the left of A
228 *
229 CALL ZGEMV( 'C', IXFRM, N, CONE, A( KBEG, 1 ), LDA,
230 $ X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 )
231 CALL ZGERC( IXFRM, N, -DCMPLX( FACTOR ), X( KBEG ), 1,
232 $ X( 2*NXFRM+1 ), 1, A( KBEG, 1 ), LDA )
233 *
234 END IF
235 *
236 IF( ITYPE.GE.2 .AND. ITYPE.LE.4 ) THEN
237 *
238 * Apply H(k)* (or H(k)') on the right of A
239 *
240 IF( ITYPE.EQ.4 ) THEN
241 CALL ZLACGV( IXFRM, X( KBEG ), 1 )
242 END IF
243 *
244 CALL ZGEMV( 'N', M, IXFRM, CONE, A( 1, KBEG ), LDA,
245 $ X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 )
246 CALL ZGERC( M, IXFRM, -DCMPLX( FACTOR ), X( 2*NXFRM+1 ), 1,
247 $ X( KBEG ), 1, A( 1, KBEG ), LDA )
248 *
249 END IF
250 30 CONTINUE
251 *
252 X( 1 ) = ZLARND( 3, ISEED )
253 XABS = ABS( X( 1 ) )
254 IF( XABS.NE.ZERO ) THEN
255 CSIGN = X( 1 ) / XABS
256 ELSE
257 CSIGN = CONE
258 END IF
259 X( 2*NXFRM ) = CSIGN
260 *
261 * Scale the matrix A by D.
262 *
263 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN
264 DO 40 IROW = 1, M
265 CALL ZSCAL( N, DCONJG( X( NXFRM+IROW ) ), A( IROW, 1 ),
266 $ LDA )
267 40 CONTINUE
268 END IF
269 *
270 IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
271 DO 50 JCOL = 1, N
272 CALL ZSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
273 50 CONTINUE
274 END IF
275 *
276 IF( ITYPE.EQ.4 ) THEN
277 DO 60 JCOL = 1, N
278 CALL ZSCAL( M, DCONJG( X( NXFRM+JCOL ) ), A( 1, JCOL ), 1 )
279 60 CONTINUE
280 END IF
281 RETURN
282 *
283 * End of ZLAROR
284 *
285 END