LAPACK  3.10.1 LAPACK: Linear Algebra PACKage
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)
XERBLA
Definition: xerbla.f:60
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERC
Definition: cgerc.f:130
subroutine claror(SIDE, INIT, M, N, A, LDA, ISEED, X, INFO)
CLAROR
Definition: claror.f:158
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