LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slaror.f
Go to the documentation of this file.
1 *> \brief \b SLAROR
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 SLAROR( 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 * REAL A( LDA, * ), X( * )
20 * ..
21 *
22 *
23 *> \par Purpose:
24 * =============
25 *>
26 *> \verbatim
27 *>
28 *> SLAROR 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 REAL 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 SLAROR to continue the same random number
113 *> sequence.
114 *> \endverbatim
115 *>
116 *> \param[out] X
117 *> \verbatim
118 *> X is REAL 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 SLARND 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 *> \ingroup real_matgen
143 *
144 * =====================================================================
145  SUBROUTINE slaror( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
146 *
147 * -- LAPACK auxiliary routine --
148 * -- LAPACK is a software package provided by Univ. of Tennessee, --
149 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150 *
151 * .. Scalar Arguments ..
152  CHARACTER INIT, SIDE
153  INTEGER INFO, LDA, M, N
154 * ..
155 * .. Array Arguments ..
156  INTEGER ISEED( 4 )
157  REAL A( LDA, * ), X( * )
158 * ..
159 *
160 * =====================================================================
161 *
162 * .. Parameters ..
163  REAL ZERO, ONE, TOOSML
164  parameter( zero = 0.0e+0, one = 1.0e+0,
165  $ toosml = 1.0e-20 )
166 * ..
167 * .. Local Scalars ..
168  INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
169  REAL FACTOR, XNORM, XNORMS
170 * ..
171 * .. External Functions ..
172  LOGICAL LSAME
173  REAL SLARND, SNRM2
174  EXTERNAL lsame, slarnd, snrm2
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL sgemv, sger, slaset, sscal, xerbla
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC abs, sign
181 * ..
182 * .. Executable Statements ..
183 *
184  info = 0
185  IF( n.EQ.0 .OR. m.EQ.0 )
186  $ RETURN
187 *
188  itype = 0
189  IF( lsame( side, 'L' ) ) THEN
190  itype = 1
191  ELSE IF( lsame( side, 'R' ) ) THEN
192  itype = 2
193  ELSE IF( lsame( side, 'C' ) .OR. lsame( side, 'T' ) ) THEN
194  itype = 3
195  END IF
196 *
197 * Check for argument errors.
198 *
199  IF( itype.EQ.0 ) THEN
200  info = -1
201  ELSE IF( m.LT.0 ) THEN
202  info = -3
203  ELSE IF( n.LT.0 .OR. ( itype.EQ.3 .AND. n.NE.m ) ) THEN
204  info = -4
205  ELSE IF( lda.LT.m ) THEN
206  info = -6
207  END IF
208  IF( info.NE.0 ) THEN
209  CALL xerbla( 'SLAROR', -info )
210  RETURN
211  END IF
212 *
213  IF( itype.EQ.1 ) THEN
214  nxfrm = m
215  ELSE
216  nxfrm = n
217  END IF
218 *
219 * Initialize A to the identity matrix if desired
220 *
221  IF( lsame( init, 'I' ) )
222  $ CALL slaset( 'Full', m, n, zero, one, a, lda )
223 *
224 * If no rotation possible, multiply by random +/-1
225 *
226 * Compute rotation by computing Householder transformations
227 * H(2), H(3), ..., H(nhouse)
228 *
229  DO 10 j = 1, nxfrm
230  x( j ) = zero
231  10 CONTINUE
232 *
233  DO 30 ixfrm = 2, nxfrm
234  kbeg = nxfrm - ixfrm + 1
235 *
236 * Generate independent normal( 0, 1 ) random numbers
237 *
238  DO 20 j = kbeg, nxfrm
239  x( j ) = slarnd( 3, iseed )
240  20 CONTINUE
241 *
242 * Generate a Householder transformation from the random vector X
243 *
244  xnorm = snrm2( ixfrm, x( kbeg ), 1 )
245  xnorms = sign( xnorm, x( kbeg ) )
246  x( kbeg+nxfrm ) = sign( one, -x( kbeg ) )
247  factor = xnorms*( xnorms+x( kbeg ) )
248  IF( abs( factor ).LT.toosml ) THEN
249  info = 1
250  CALL xerbla( 'SLAROR', info )
251  RETURN
252  ELSE
253  factor = one / factor
254  END IF
255  x( kbeg ) = x( kbeg ) + xnorms
256 *
257 * Apply Householder transformation to A
258 *
259  IF( itype.EQ.1 .OR. itype.EQ.3 ) THEN
260 *
261 * Apply H(k) from the left.
262 *
263  CALL sgemv( 'T', ixfrm, n, one, a( kbeg, 1 ), lda,
264  $ x( kbeg ), 1, zero, x( 2*nxfrm+1 ), 1 )
265  CALL sger( ixfrm, n, -factor, x( kbeg ), 1, x( 2*nxfrm+1 ),
266  $ 1, a( kbeg, 1 ), lda )
267 *
268  END IF
269 *
270  IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
271 *
272 * Apply H(k) from the right.
273 *
274  CALL sgemv( 'N', m, ixfrm, one, a( 1, kbeg ), lda,
275  $ x( kbeg ), 1, zero, x( 2*nxfrm+1 ), 1 )
276  CALL sger( m, ixfrm, -factor, x( 2*nxfrm+1 ), 1, x( kbeg ),
277  $ 1, a( 1, kbeg ), lda )
278 *
279  END IF
280  30 CONTINUE
281 *
282  x( 2*nxfrm ) = sign( one, slarnd( 3, iseed ) )
283 *
284 * Scale the matrix A by D.
285 *
286  IF( itype.EQ.1 .OR. itype.EQ.3 ) THEN
287  DO 40 irow = 1, m
288  CALL sscal( n, x( nxfrm+irow ), a( irow, 1 ), lda )
289  40 CONTINUE
290  END IF
291 *
292  IF( itype.EQ.2 .OR. itype.EQ.3 ) THEN
293  DO 50 jcol = 1, n
294  CALL sscal( m, x( nxfrm+jcol ), a( 1, jcol ), 1 )
295  50 CONTINUE
296  END IF
297  RETURN
298 *
299 * End of SLAROR
300 *
301  END
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slaror(SIDE, INIT, M, N, A, LDA, ISEED, X, INFO)
SLAROR
Definition: slaror.f:146
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:130
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156