LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
claror.f
Go to the documentation of this file.
1*> \brief \b CLAROR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE CLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
12*
13* .. Scalar Arguments ..
14* CHARACTER INIT, SIDE
15* INTEGER INFO, LDA, M, N
16* ..
17* .. Array Arguments ..
18* INTEGER ISEED( 4 )
19* COMPLEX A( LDA, * ), X( * )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> CLAROR pre- or post-multiplies an M by N matrix A by a random
29*> unitary matrix U, overwriting A. A may optionally be
30*> initialized to the identity matrix before multiplying by U.
31*> U is generated using the method of G.W. Stewart
32*> ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ).
33*> (BLAS-2 version)
34*> \endverbatim
35*
36* Arguments:
37* ==========
38*
39*> \param[in] SIDE
40*> \verbatim
41*> SIDE is CHARACTER*1
42*> SIDE specifies whether A is multiplied on the left or right
43*> by U.
44*> SIDE = 'L' Multiply A on the left (premultiply) by U
45*> SIDE = 'R' Multiply A on the right (postmultiply) by UC> SIDE = 'C' Multiply A on the left by U and the right by UC> SIDE = 'T' Multiply A on the left by U and the right by U'
46*> Not modified.
47*> \endverbatim
48*>
49*> \param[in] INIT
50*> \verbatim
51*> INIT is CHARACTER*1
52*> INIT specifies whether or not A should be initialized to
53*> the identity matrix.
54*> INIT = 'I' Initialize A to (a section of) the
55*> identity matrix before applying U.
56*> INIT = 'N' No initialization. Apply U to the
57*> input matrix A.
58*>
59*> INIT = 'I' may be used to generate square (i.e., unitary)
60*> or rectangular orthogonal matrices (orthogonality being
61*> in the sense of CDOTC):
62*>
63*> For square matrices, M=N, and SIDE many be either 'L' or
64*> 'R'; the rows will be orthogonal to each other, as will the
65*> columns.
66*> For rectangular matrices where M < N, SIDE = 'R' will
67*> produce a dense matrix whose rows will be orthogonal and
68*> whose columns will not, while SIDE = 'L' will produce a
69*> matrix whose rows will be orthogonal, and whose first M
70*> columns will be orthogonal, the remaining columns being
71*> zero.
72*> For matrices where M > N, just use the previous
73*> explanation, interchanging 'L' and 'R' and "rows" and
74*> "columns".
75*>
76*> Not modified.
77*> \endverbatim
78*>
79*> \param[in] M
80*> \verbatim
81*> M is INTEGER
82*> Number of rows of A. Not modified.
83*> \endverbatim
84*>
85*> \param[in] N
86*> \verbatim
87*> N is INTEGER
88*> Number of columns of A. Not modified.
89*> \endverbatim
90*>
91*> \param[in,out] A
92*> \verbatim
93*> A is COMPLEX array, dimension ( LDA, N )
94*> Input and output array. Overwritten by U A ( if SIDE = 'L' )
95*> or by A U ( if SIDE = 'R' )
96*> or by U A U* ( if SIDE = 'C')
97*> or by U A U' ( if SIDE = 'T') on exit.
98*> \endverbatim
99*>
100*> \param[in] LDA
101*> \verbatim
102*> LDA is INTEGER
103*> Leading dimension of A. Must be at least MAX ( 1, M ).
104*> Not modified.
105*> \endverbatim
106*>
107*> \param[in,out] ISEED
108*> \verbatim
109*> ISEED is INTEGER array, dimension ( 4 )
110*> On entry ISEED specifies the seed of the random number
111*> generator. The array elements should be between 0 and 4095;
112*> if not they will be reduced mod 4096. Also, ISEED(4) must
113*> be odd. The random number generator uses a linear
114*> congruential sequence limited to small integers, and so
115*> should produce machine independent random numbers. The
116*> values of ISEED are changed on exit, and can be used in the
117*> next call to CLAROR to continue the same random number
118*> sequence.
119*> Modified.
120*> \endverbatim
121*>
122*> \param[out] X
123*> \verbatim
124*> X is COMPLEX array, dimension ( 3*MAX( M, N ) )
125*> Workspace. Of length:
126*> 2*M + N if SIDE = 'L',
127*> 2*N + M if SIDE = 'R',
128*> 3*N if SIDE = 'C' or 'T'.
129*> Modified.
130*> \endverbatim
131*>
132*> \param[out] INFO
133*> \verbatim
134*> INFO is INTEGER
135*> An error flag. It is set to:
136*> 0 if no error.
137*> 1 if CLARND returned a bad random number (installation
138*> problem)
139*> -1 if SIDE is not L, R, C, or T.
140*> -3 if M is negative.
141*> -4 if N is negative or if SIDE is C or T and N is not equal
142*> to M.
143*> -6 if LDA is less than M.
144*> \endverbatim
145*
146* Authors:
147* ========
148*
149*> \author Univ. of Tennessee
150*> \author Univ. of California Berkeley
151*> \author Univ. of Colorado Denver
152*> \author NAG Ltd.
153*
154*> \ingroup complex_matgen
155*
156* =====================================================================
157 SUBROUTINE claror( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
158*
159* -- LAPACK auxiliary routine --
160* -- LAPACK is a software package provided by Univ. of Tennessee, --
161* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162*
163* .. Scalar Arguments ..
164 CHARACTER INIT, SIDE
165 INTEGER INFO, LDA, M, N
166* ..
167* .. Array Arguments ..
168 INTEGER ISEED( 4 )
169 COMPLEX A( LDA, * ), X( * )
170* ..
171*
172* =====================================================================
173*
174* .. Parameters ..
175 REAL ZERO, ONE, TOOSML
176 parameter( zero = 0.0e+0, one = 1.0e+0,
177 $ toosml = 1.0e-20 )
178 COMPLEX CZERO, CONE
179 parameter( czero = ( 0.0e+0, 0.0e+0 ),
180 $ cone = ( 1.0e+0, 0.0e+0 ) )
181* ..
182* .. Local Scalars ..
183 INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
184 REAL FACTOR, XABS, XNORM
185 COMPLEX CSIGN, XNORMS
186* ..
187* .. External Functions ..
188 LOGICAL LSAME
189 REAL SCNRM2
190 COMPLEX CLARND
191 EXTERNAL lsame, scnrm2, clarnd
192* ..
193* .. External Subroutines ..
194 EXTERNAL cgemv, cgerc, clacgv, claset, cscal, xerbla
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC abs, cmplx, conjg
198* ..
199* .. Executable Statements ..
200*
201 info = 0
202 IF( n.EQ.0 .OR. m.EQ.0 )
203 $ RETURN
204*
205 itype = 0
206 IF( lsame( side, 'L' ) ) THEN
207 itype = 1
208 ELSE IF( lsame( side, 'R' ) ) THEN
209 itype = 2
210 ELSE IF( lsame( side, 'C' ) ) THEN
211 itype = 3
212 ELSE IF( lsame( side, 'T' ) ) THEN
213 itype = 4
214 END IF
215*
216* Check for argument errors.
217*
218 IF( itype.EQ.0 ) THEN
219 info = -1
220 ELSE IF( m.LT.0 ) THEN
221 info = -3
222 ELSE IF( n.LT.0 .OR. ( itype.EQ.3 .AND. n.NE.m ) ) THEN
223 info = -4
224 ELSE IF( lda.LT.m ) THEN
225 info = -6
226 END IF
227 IF( info.NE.0 ) THEN
228 CALL xerbla( 'CLAROR', -info )
229 RETURN
230 END IF
231*
232 IF( itype.EQ.1 ) THEN
233 nxfrm = m
234 ELSE
235 nxfrm = n
236 END IF
237*
238* Initialize A to the identity matrix if desired
239*
240 IF( lsame( init, 'I' ) )
241 $ CALL claset( 'Full', m, n, czero, cone, a, lda )
242*
243* If no rotation possible, still multiply by
244* a random complex number from the circle |x| = 1
245*
246* 2) Compute Rotation by computing Householder
247* Transformations H(2), H(3), ..., H(n). Note that the
248* order in which they are computed is irrelevant.
249*
250 DO 40 j = 1, nxfrm
251 x( j ) = czero
252 40 CONTINUE
253*
254 DO 60 ixfrm = 2, nxfrm
255 kbeg = nxfrm - ixfrm + 1
256*
257* Generate independent normal( 0, 1 ) random numbers
258*
259 DO 50 j = kbeg, nxfrm
260 x( j ) = clarnd( 3, iseed )
261 50 CONTINUE
262*
263* Generate a Householder transformation from the random vector X
264*
265 xnorm = scnrm2( ixfrm, x( kbeg ), 1 )
266 xabs = abs( x( kbeg ) )
267 IF( xabs.NE.czero ) THEN
268 csign = x( kbeg ) / xabs
269 ELSE
270 csign = cone
271 END IF
272 xnorms = csign*xnorm
273 x( nxfrm+kbeg ) = -csign
274 factor = xnorm*( xnorm+xabs )
275 IF( abs( factor ).LT.toosml ) THEN
276 info = 1
277 CALL xerbla( 'CLAROR', -info )
278 RETURN
279 ELSE
280 factor = one / factor
281 END IF
282 x( kbeg ) = x( kbeg ) + xnorms
283*
284* Apply Householder transformation to A
285*
286 IF( itype.EQ.1 .OR. itype.EQ.3 .OR. itype.EQ.4 ) THEN
287*
288* Apply H(k) on the left of A
289*
290 CALL cgemv( 'C', ixfrm, n, cone, a( kbeg, 1 ), lda,
291 $ x( kbeg ), 1, czero, x( 2*nxfrm+1 ), 1 )
292 CALL cgerc( ixfrm, n, -cmplx( factor ), x( kbeg ), 1,
293 $ x( 2*nxfrm+1 ), 1, a( kbeg, 1 ), lda )
294*
295 END IF
296*
297 IF( itype.GE.2 .AND. itype.LE.4 ) THEN
298*
299* Apply H(k)* (or H(k)') on the right of A
300*
301 IF( itype.EQ.4 ) THEN
302 CALL clacgv( ixfrm, x( kbeg ), 1 )
303 END IF
304*
305 CALL cgemv( 'N', m, ixfrm, cone, a( 1, kbeg ), lda,
306 $ x( kbeg ), 1, czero, x( 2*nxfrm+1 ), 1 )
307 CALL cgerc( m, ixfrm, -cmplx( factor ), x( 2*nxfrm+1 ), 1,
308 $ x( kbeg ), 1, a( 1, kbeg ), lda )
309*
310 END IF
311 60 CONTINUE
312*
313 x( 1 ) = clarnd( 3, iseed )
314 xabs = abs( x( 1 ) )
315 IF( xabs.NE.zero ) THEN
316 csign = x( 1 ) / xabs
317 ELSE
318 csign = cone
319 END IF
320 x( 2*nxfrm ) = csign
321*
322* Scale the matrix A by D.
323*
324 IF( itype.EQ.1 .OR. itype.EQ.3 .OR. itype.EQ.4 ) THEN
325 DO 70 irow = 1, m
326 CALL cscal( n, conjg( x( nxfrm+irow ) ), a( irow, 1 ), lda )
327 70 CONTINUE
328 END IF
329*
330 IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
331 DO 80 jcol = 1, n
332 CALL cscal( m, x( nxfrm+jcol ), a( 1, jcol ), 1 )
333 80 CONTINUE
334 END IF
335*
336 IF( itype.EQ.4 ) THEN
337 DO 90 jcol = 1, n
338 CALL cscal( m, conjg( x( nxfrm+jcol ) ), a( 1, jcol ), 1 )
339 90 CONTINUE
340 END IF
341 RETURN
342*
343* End of CLAROR
344*
345 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine claror(side, init, m, n, a, lda, iseed, x, info)
CLAROR
Definition claror.f:158
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine cgerc(m, n, alpha, x, incx, y, incy, a, lda)
CGERC
Definition cgerc.f:130
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78