LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlaror.f
Go to the documentation of this file.
1 *> \brief \b ZLAROR
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 ZLAROR( 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*16 A( LDA, * ), X( * )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> ZLAROR 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 ZDOTC):
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*16 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 ZLAROR to continue the same random number
118 *> sequence.
119 *> Modified.
120 *> \endverbatim
121 *>
122 *> \param[out] X
123 *> \verbatim
124 *> X is COMPLEX*16 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 ZLARND 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 *> \date November 2011
155 *
156 *> \ingroup complex16_matgen
157 *
158 * =====================================================================
159  SUBROUTINE zlaror( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
160 *
161 * -- LAPACK auxiliary routine (version 3.4.0) --
162 * -- LAPACK is a software package provided by Univ. of Tennessee, --
163 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164 * November 2011
165 *
166 * .. Scalar Arguments ..
167  CHARACTER init, side
168  INTEGER info, lda, m, n
169 * ..
170 * .. Array Arguments ..
171  INTEGER iseed( 4 )
172  COMPLEX*16 a( lda, * ), x( * )
173 * ..
174 *
175 * =====================================================================
176 *
177 * .. Parameters ..
178  DOUBLE PRECISION zero, one, toosml
179  parameter( zero = 0.0d+0, one = 1.0d+0,
180  $ toosml = 1.0d-20 )
181  COMPLEX*16 czero, cone
182  parameter( czero = ( 0.0d+0, 0.0d+0 ),
183  $ cone = ( 1.0d+0, 0.0d+0 ) )
184 * ..
185 * .. Local Scalars ..
186  INTEGER irow, itype, ixfrm, j, jcol, kbeg, nxfrm
187  DOUBLE PRECISION factor, xabs, xnorm
188  COMPLEX*16 csign, xnorms
189 * ..
190 * .. External Functions ..
191  LOGICAL lsame
192  DOUBLE PRECISION dznrm2
193  COMPLEX*16 zlarnd
194  EXTERNAL lsame, dznrm2, zlarnd
195 * ..
196 * .. External Subroutines ..
197  EXTERNAL xerbla, zgemv, zgerc, zlacgv, zlaset, zscal
198 * ..
199 * .. Intrinsic Functions ..
200  INTRINSIC abs, dcmplx, dconjg
201 * ..
202 * .. Executable Statements ..
203 *
204  info = 0
205  IF( n.EQ.0 .OR. m.EQ.0 )
206  $ return
207 *
208  itype = 0
209  IF( lsame( side, 'L' ) ) THEN
210  itype = 1
211  ELSE IF( lsame( side, 'R' ) ) THEN
212  itype = 2
213  ELSE IF( lsame( side, 'C' ) ) THEN
214  itype = 3
215  ELSE IF( lsame( side, 'T' ) ) THEN
216  itype = 4
217  END IF
218 *
219 * Check for argument errors.
220 *
221  IF( itype.EQ.0 ) THEN
222  info = -1
223  ELSE IF( m.LT.0 ) THEN
224  info = -3
225  ELSE IF( n.LT.0 .OR. ( itype.EQ.3 .AND. n.NE.m ) ) THEN
226  info = -4
227  ELSE IF( lda.LT.m ) THEN
228  info = -6
229  END IF
230  IF( info.NE.0 ) THEN
231  CALL xerbla( 'ZLAROR', -info )
232  return
233  END IF
234 *
235  IF( itype.EQ.1 ) THEN
236  nxfrm = m
237  ELSE
238  nxfrm = n
239  END IF
240 *
241 * Initialize A to the identity matrix if desired
242 *
243  IF( lsame( init, 'I' ) )
244  $ CALL zlaset( 'Full', m, n, czero, cone, a, lda )
245 *
246 * If no rotation possible, still multiply by
247 * a random complex number from the circle |x| = 1
248 *
249 * 2) Compute Rotation by computing Householder
250 * Transformations H(2), H(3), ..., H(n). Note that the
251 * order in which they are computed is irrelevant.
252 *
253  DO 10 j = 1, nxfrm
254  x( j ) = czero
255  10 continue
256 *
257  DO 30 ixfrm = 2, nxfrm
258  kbeg = nxfrm - ixfrm + 1
259 *
260 * Generate independent normal( 0, 1 ) random numbers
261 *
262  DO 20 j = kbeg, nxfrm
263  x( j ) = zlarnd( 3, iseed )
264  20 continue
265 *
266 * Generate a Householder transformation from the random vector X
267 *
268  xnorm = dznrm2( ixfrm, x( kbeg ), 1 )
269  xabs = abs( x( kbeg ) )
270  IF( xabs.NE.czero ) THEN
271  csign = x( kbeg ) / xabs
272  ELSE
273  csign = cone
274  END IF
275  xnorms = csign*xnorm
276  x( nxfrm+kbeg ) = -csign
277  factor = xnorm*( xnorm+xabs )
278  IF( abs( factor ).LT.toosml ) THEN
279  info = 1
280  CALL xerbla( 'ZLAROR', -info )
281  return
282  ELSE
283  factor = one / factor
284  END IF
285  x( kbeg ) = x( kbeg ) + xnorms
286 *
287 * Apply Householder transformation to A
288 *
289  IF( itype.EQ.1 .OR. itype.EQ.3 .OR. itype.EQ.4 ) THEN
290 *
291 * Apply H(k) on the left of A
292 *
293  CALL zgemv( 'C', ixfrm, n, cone, a( kbeg, 1 ), lda,
294  $ x( kbeg ), 1, czero, x( 2*nxfrm+1 ), 1 )
295  CALL zgerc( ixfrm, n, -dcmplx( factor ), x( kbeg ), 1,
296  $ x( 2*nxfrm+1 ), 1, a( kbeg, 1 ), lda )
297 *
298  END IF
299 *
300  IF( itype.GE.2 .AND. itype.LE.4 ) THEN
301 *
302 * Apply H(k)* (or H(k)') on the right of A
303 *
304  IF( itype.EQ.4 ) THEN
305  CALL zlacgv( ixfrm, x( kbeg ), 1 )
306  END IF
307 *
308  CALL zgemv( 'N', m, ixfrm, cone, a( 1, kbeg ), lda,
309  $ x( kbeg ), 1, czero, x( 2*nxfrm+1 ), 1 )
310  CALL zgerc( m, ixfrm, -dcmplx( factor ), x( 2*nxfrm+1 ), 1,
311  $ x( kbeg ), 1, a( 1, kbeg ), lda )
312 *
313  END IF
314  30 continue
315 *
316  x( 1 ) = zlarnd( 3, iseed )
317  xabs = abs( x( 1 ) )
318  IF( xabs.NE.zero ) THEN
319  csign = x( 1 ) / xabs
320  ELSE
321  csign = cone
322  END IF
323  x( 2*nxfrm ) = csign
324 *
325 * Scale the matrix A by D.
326 *
327  IF( itype.EQ.1 .OR. itype.EQ.3 .OR. itype.EQ.4 ) THEN
328  DO 40 irow = 1, m
329  CALL zscal( n, dconjg( x( nxfrm+irow ) ), a( irow, 1 ),
330  $ lda )
331  40 continue
332  END IF
333 *
334  IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
335  DO 50 jcol = 1, n
336  CALL zscal( m, x( nxfrm+jcol ), a( 1, jcol ), 1 )
337  50 continue
338  END IF
339 *
340  IF( itype.EQ.4 ) THEN
341  DO 60 jcol = 1, n
342  CALL zscal( m, dconjg( x( nxfrm+jcol ) ), a( 1, jcol ), 1 )
343  60 continue
344  END IF
345  return
346 *
347 * End of ZLAROR
348 *
349  END