LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
dlaror.f
Go to the documentation of this file.
1 *> \brief \b DLAROR
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 DLAROR( 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 * DOUBLE PRECISION A( LDA, * ), X( * )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> DLAROR pre- or post-multiplies an M by N matrix A by a random
29 *> orthogonal matrix U, overwriting A. A may optionally be initialized
30 *> to the identity matrix before multiplying by U. U is generated using
31 *> the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
32 *> \endverbatim
33 *
34 * Arguments:
35 * ==========
36 *
37 *> \param[in] SIDE
38 *> \verbatim
39 *> SIDE is CHARACTER*1
40 *> Specifies whether A is multiplied on the left or right by U.
41 *> = 'L': Multiply A on the left (premultiply) by U
42 *> = 'R': Multiply A on the right (postmultiply) by U'
43 *> = 'C' or 'T': Multiply A on the left by U and the right
44 *> by U' (Here, U' means U-transpose.)
45 *> \endverbatim
46 *>
47 *> \param[in] INIT
48 *> \verbatim
49 *> INIT is CHARACTER*1
50 *> Specifies whether or not A should be initialized to the
51 *> identity matrix.
52 *> = 'I': Initialize A to (a section of) the identity matrix
53 *> before applying U.
54 *> = 'N': No initialization. Apply U to the input matrix A.
55 *>
56 *> INIT = 'I' may be used to generate square or rectangular
57 *> orthogonal matrices:
58 *>
59 *> For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
60 *> to each other, as will the columns.
61 *>
62 *> If M < N, SIDE = 'R' produces a dense matrix whose rows are
63 *> orthogonal and whose columns are not, while SIDE = 'L'
64 *> produces a matrix whose rows are orthogonal, and whose first
65 *> M columns are orthogonal, and whose remaining columns are
66 *> zero.
67 *>
68 *> If M > N, SIDE = 'L' produces a dense matrix whose columns
69 *> are orthogonal and whose rows are not, while SIDE = 'R'
70 *> produces a matrix whose columns are orthogonal, and whose
71 *> first M rows are orthogonal, and whose remaining rows are
72 *> zero.
73 *> \endverbatim
74 *>
75 *> \param[in] M
76 *> \verbatim
77 *> M is INTEGER
78 *> The number of rows of A.
79 *> \endverbatim
80 *>
81 *> \param[in] N
82 *> \verbatim
83 *> N is INTEGER
84 *> The number of columns of A.
85 *> \endverbatim
86 *>
87 *> \param[in,out] A
88 *> \verbatim
89 *> A is DOUBLE PRECISION array, dimension (LDA, N)
90 *> On entry, the array A.
91 *> On exit, overwritten by U A ( if SIDE = 'L' ),
92 *> or by A U ( if SIDE = 'R' ),
93 *> or by U A U' ( if SIDE = 'C' or 'T').
94 *> \endverbatim
95 *>
96 *> \param[in] LDA
97 *> \verbatim
98 *> LDA is INTEGER
99 *> The leading dimension of the array A. LDA >= max(1,M).
100 *> \endverbatim
101 *>
102 *> \param[in,out] ISEED
103 *> \verbatim
104 *> ISEED is INTEGER array, dimension (4)
105 *> On entry ISEED specifies the seed of the random number
106 *> generator. The array elements should be between 0 and 4095;
107 *> if not they will be reduced mod 4096. Also, ISEED(4) must
108 *> be odd. The random number generator uses a linear
109 *> congruential sequence limited to small integers, and so
110 *> should produce machine independent random numbers. The
111 *> values of ISEED are changed on exit, and can be used in the
112 *> next call to DLAROR to continue the same random number
113 *> sequence.
114 *> \endverbatim
115 *>
116 *> \param[out] X
117 *> \verbatim
118 *> X is DOUBLE PRECISION array, dimension (3*MAX( M, N ))
119 *> Workspace of length
120 *> 2*M + N if SIDE = 'L',
121 *> 2*N + M if SIDE = 'R',
122 *> 3*N if SIDE = 'C' or 'T'.
123 *> \endverbatim
124 *>
125 *> \param[out] INFO
126 *> \verbatim
127 *> INFO is INTEGER
128 *> An error flag. It is set to:
129 *> = 0: normal return
130 *> < 0: if INFO = -k, the k-th argument had an illegal value
131 *> = 1: if the random numbers generated by DLARND are bad.
132 *> \endverbatim
133 *
134 * Authors:
135 * ========
136 *
137 *> \author Univ. of Tennessee
138 *> \author Univ. of California Berkeley
139 *> \author Univ. of Colorado Denver
140 *> \author NAG Ltd.
141 *
142 *> \date November 2011
143 *
144 *> \ingroup double_matgen
145 *
146 * =====================================================================
147  SUBROUTINE dlaror( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
148 *
149 * -- LAPACK auxiliary routine (version 3.4.0) --
150 * -- LAPACK is a software package provided by Univ. of Tennessee, --
151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152 * November 2011
153 *
154 * .. Scalar Arguments ..
155  CHARACTER init, side
156  INTEGER info, lda, m, n
157 * ..
158 * .. Array Arguments ..
159  INTEGER iseed( 4 )
160  DOUBLE PRECISION a( lda, * ), x( * )
161 * ..
162 *
163 * =====================================================================
164 *
165 * .. Parameters ..
166  DOUBLE PRECISION zero, one, toosml
167  parameter( zero = 0.0d+0, one = 1.0d+0,
168  \$ toosml = 1.0d-20 )
169 * ..
170 * .. Local Scalars ..
171  INTEGER irow, itype, ixfrm, j, jcol, kbeg, nxfrm
172  DOUBLE PRECISION factor, xnorm, xnorms
173 * ..
174 * .. External Functions ..
175  LOGICAL lsame
176  DOUBLE PRECISION dlarnd, dnrm2
177  EXTERNAL lsame, dlarnd, dnrm2
178 * ..
179 * .. External Subroutines ..
180  EXTERNAL dgemv, dger, dlaset, dscal, xerbla
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC abs, sign
184 * ..
185 * .. Executable Statements ..
186 *
187  info = 0
188  IF( n.EQ.0 .OR. m.EQ.0 )
189  \$ return
190 *
191  itype = 0
192  IF( lsame( side, 'L' ) ) THEN
193  itype = 1
194  ELSE IF( lsame( side, 'R' ) ) THEN
195  itype = 2
196  ELSE IF( lsame( side, 'C' ) .OR. lsame( side, 'T' ) ) THEN
197  itype = 3
198  END IF
199 *
200 * Check for argument errors.
201 *
202  IF( itype.EQ.0 ) THEN
203  info = -1
204  ELSE IF( m.LT.0 ) THEN
205  info = -3
206  ELSE IF( n.LT.0 .OR. ( itype.EQ.3 .AND. n.NE.m ) ) THEN
207  info = -4
208  ELSE IF( lda.LT.m ) THEN
209  info = -6
210  END IF
211  IF( info.NE.0 ) THEN
212  CALL xerbla( 'DLAROR', -info )
213  return
214  END IF
215 *
216  IF( itype.EQ.1 ) THEN
217  nxfrm = m
218  ELSE
219  nxfrm = n
220  END IF
221 *
222 * Initialize A to the identity matrix if desired
223 *
224  IF( lsame( init, 'I' ) )
225  \$ CALL dlaset( 'Full', m, n, zero, one, a, lda )
226 *
227 * If no rotation possible, multiply by random +/-1
228 *
229 * Compute rotation by computing Householder transformations
230 * H(2), H(3), ..., H(nhouse)
231 *
232  DO 10 j = 1, nxfrm
233  x( j ) = zero
234  10 continue
235 *
236  DO 30 ixfrm = 2, nxfrm
237  kbeg = nxfrm - ixfrm + 1
238 *
239 * Generate independent normal( 0, 1 ) random numbers
240 *
241  DO 20 j = kbeg, nxfrm
242  x( j ) = dlarnd( 3, iseed )
243  20 continue
244 *
245 * Generate a Householder transformation from the random vector X
246 *
247  xnorm = dnrm2( ixfrm, x( kbeg ), 1 )
248  xnorms = sign( xnorm, x( kbeg ) )
249  x( kbeg+nxfrm ) = sign( one, -x( kbeg ) )
250  factor = xnorms*( xnorms+x( kbeg ) )
251  IF( abs( factor ).LT.toosml ) THEN
252  info = 1
253  CALL xerbla( 'DLAROR', info )
254  return
255  ELSE
256  factor = one / factor
257  END IF
258  x( kbeg ) = x( kbeg ) + xnorms
259 *
260 * Apply Householder transformation to A
261 *
262  IF( itype.EQ.1 .OR. itype.EQ.3 ) THEN
263 *
264 * Apply H(k) from the left.
265 *
266  CALL dgemv( 'T', ixfrm, n, one, a( kbeg, 1 ), lda,
267  \$ x( kbeg ), 1, zero, x( 2*nxfrm+1 ), 1 )
268  CALL dger( ixfrm, n, -factor, x( kbeg ), 1, x( 2*nxfrm+1 ),
269  \$ 1, a( kbeg, 1 ), lda )
270 *
271  END IF
272 *
273  IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
274 *
275 * Apply H(k) from the right.
276 *
277  CALL dgemv( 'N', m, ixfrm, one, a( 1, kbeg ), lda,
278  \$ x( kbeg ), 1, zero, x( 2*nxfrm+1 ), 1 )
279  CALL dger( m, ixfrm, -factor, x( 2*nxfrm+1 ), 1, x( kbeg ),
280  \$ 1, a( 1, kbeg ), lda )
281 *
282  END IF
283  30 continue
284 *
285  x( 2*nxfrm ) = sign( one, dlarnd( 3, iseed ) )
286 *
287 * Scale the matrix A by D.
288 *
289  IF( itype.EQ.1 .OR. itype.EQ.3 ) THEN
290  DO 40 irow = 1, m
291  CALL dscal( n, x( nxfrm+irow ), a( irow, 1 ), lda )
292  40 continue
293  END IF
294 *
295  IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
296  DO 50 jcol = 1, n
297  CALL dscal( m, x( nxfrm+jcol ), a( 1, jcol ), 1 )
298  50 continue
299  END IF
300  return
301 *
302 * End of DLAROR
303 *
304  END