LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlagge.f
Go to the documentation of this file.
1*> \brief \b DLAGGE
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 DLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO )
12*
13* .. Scalar Arguments ..
14* INTEGER INFO, KL, KU, LDA, M, N
15* ..
16* .. Array Arguments ..
17* INTEGER ISEED( 4 )
18* DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
19* ..
20*
21*
22*> \par Purpose:
23* =============
24*>
25*> \verbatim
26*>
27*> DLAGGE generates a real general m by n matrix A, by pre- and post-
28*> multiplying a real diagonal matrix D with random orthogonal matrices:
29*> A = U*D*V. The lower and upper bandwidths may then be reduced to
30*> kl and ku by additional orthogonal transformations.
31*> \endverbatim
32*
33* Arguments:
34* ==========
35*
36*> \param[in] M
37*> \verbatim
38*> M is INTEGER
39*> The number of rows of the matrix A. M >= 0.
40*> \endverbatim
41*>
42*> \param[in] N
43*> \verbatim
44*> N is INTEGER
45*> The number of columns of the matrix A. N >= 0.
46*> \endverbatim
47*>
48*> \param[in] KL
49*> \verbatim
50*> KL is INTEGER
51*> The number of nonzero subdiagonals within the band of A.
52*> 0 <= KL <= M-1.
53*> \endverbatim
54*>
55*> \param[in] KU
56*> \verbatim
57*> KU is INTEGER
58*> The number of nonzero superdiagonals within the band of A.
59*> 0 <= KU <= N-1.
60*> \endverbatim
61*>
62*> \param[in] D
63*> \verbatim
64*> D is DOUBLE PRECISION array, dimension (min(M,N))
65*> The diagonal elements of the diagonal matrix D.
66*> \endverbatim
67*>
68*> \param[out] A
69*> \verbatim
70*> A is DOUBLE PRECISION array, dimension (LDA,N)
71*> The generated m by n matrix A.
72*> \endverbatim
73*>
74*> \param[in] LDA
75*> \verbatim
76*> LDA is INTEGER
77*> The leading dimension of the array A. LDA >= M.
78*> \endverbatim
79*>
80*> \param[in,out] ISEED
81*> \verbatim
82*> ISEED is INTEGER array, dimension (4)
83*> On entry, the seed of the random number generator; the array
84*> elements must be between 0 and 4095, and ISEED(4) must be
85*> odd.
86*> On exit, the seed is updated.
87*> \endverbatim
88*>
89*> \param[out] WORK
90*> \verbatim
91*> WORK is DOUBLE PRECISION array, dimension (M+N)
92*> \endverbatim
93*>
94*> \param[out] INFO
95*> \verbatim
96*> INFO is INTEGER
97*> = 0: successful exit
98*> < 0: if INFO = -i, the i-th argument had an illegal value
99*> \endverbatim
100*
101* Authors:
102* ========
103*
104*> \author Univ. of Tennessee
105*> \author Univ. of California Berkeley
106*> \author Univ. of Colorado Denver
107*> \author NAG Ltd.
108*
109*> \ingroup double_matgen
110*
111* =====================================================================
112 SUBROUTINE dlagge( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO )
113*
114* -- LAPACK auxiliary routine --
115* -- LAPACK is a software package provided by Univ. of Tennessee, --
116* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
117*
118* .. Scalar Arguments ..
119 INTEGER INFO, KL, KU, LDA, M, N
120* ..
121* .. Array Arguments ..
122 INTEGER ISEED( 4 )
123 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
124* ..
125*
126* =====================================================================
127*
128* .. Parameters ..
129 DOUBLE PRECISION ZERO, ONE
130 parameter( zero = 0.0d+0, one = 1.0d+0 )
131* ..
132* .. Local Scalars ..
133 INTEGER I, J
134 DOUBLE PRECISION TAU, WA, WB, WN
135* ..
136* .. External Subroutines ..
137 EXTERNAL dgemv, dger, dlarnv, dscal, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max, min, sign
141* ..
142* .. External Functions ..
143 DOUBLE PRECISION DNRM2
144 EXTERNAL dnrm2
145* ..
146* .. Executable Statements ..
147*
148* Test the input arguments
149*
150 info = 0
151 IF( m.LT.0 ) THEN
152 info = -1
153 ELSE IF( n.LT.0 ) THEN
154 info = -2
155 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
156 info = -3
157 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
158 info = -4
159 ELSE IF( lda.LT.max( 1, m ) ) THEN
160 info = -7
161 END IF
162 IF( info.LT.0 ) THEN
163 CALL xerbla( 'DLAGGE', -info )
164 RETURN
165 END IF
166*
167* initialize A to diagonal matrix
168*
169 DO 20 j = 1, n
170 DO 10 i = 1, m
171 a( i, j ) = zero
172 10 CONTINUE
173 20 CONTINUE
174 DO 30 i = 1, min( m, n )
175 a( i, i ) = d( i )
176 30 CONTINUE
177*
178* Quick exit if the user wants a diagonal matrix
179*
180 IF(( kl .EQ. 0 ).AND.( ku .EQ. 0)) RETURN
181*
182* pre- and post-multiply A by random orthogonal matrices
183*
184 DO 40 i = min( m, n ), 1, -1
185 IF( i.LT.m ) THEN
186*
187* generate random reflection
188*
189 CALL dlarnv( 3, iseed, m-i+1, work )
190 wn = dnrm2( m-i+1, work, 1 )
191 wa = sign( wn, work( 1 ) )
192 IF( wn.EQ.zero ) THEN
193 tau = zero
194 ELSE
195 wb = work( 1 ) + wa
196 CALL dscal( m-i, one / wb, work( 2 ), 1 )
197 work( 1 ) = one
198 tau = wb / wa
199 END IF
200*
201* multiply A(i:m,i:n) by random reflection from the left
202*
203 CALL dgemv( 'Transpose', m-i+1, n-i+1, one, a( i, i ), lda,
204 $ work, 1, zero, work( m+1 ), 1 )
205 CALL dger( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
206 $ a( i, i ), lda )
207 END IF
208 IF( i.LT.n ) THEN
209*
210* generate random reflection
211*
212 CALL dlarnv( 3, iseed, n-i+1, work )
213 wn = dnrm2( n-i+1, work, 1 )
214 wa = sign( wn, work( 1 ) )
215 IF( wn.EQ.zero ) THEN
216 tau = zero
217 ELSE
218 wb = work( 1 ) + wa
219 CALL dscal( n-i, one / wb, work( 2 ), 1 )
220 work( 1 ) = one
221 tau = wb / wa
222 END IF
223*
224* multiply A(i:m,i:n) by random reflection from the right
225*
226 CALL dgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
227 $ lda, work, 1, zero, work( n+1 ), 1 )
228 CALL dger( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
229 $ a( i, i ), lda )
230 END IF
231 40 CONTINUE
232*
233* Reduce number of subdiagonals to KL and number of superdiagonals
234* to KU
235*
236 DO 70 i = 1, max( m-1-kl, n-1-ku )
237 IF( kl.LE.ku ) THEN
238*
239* annihilate subdiagonal elements first (necessary if KL = 0)
240*
241 IF( i.LE.min( m-1-kl, n ) ) THEN
242*
243* generate reflection to annihilate A(kl+i+1:m,i)
244*
245 wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
246 wa = sign( wn, a( kl+i, i ) )
247 IF( wn.EQ.zero ) THEN
248 tau = zero
249 ELSE
250 wb = a( kl+i, i ) + wa
251 CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
252 a( kl+i, i ) = one
253 tau = wb / wa
254 END IF
255*
256* apply reflection to A(kl+i:m,i+1:n) from the left
257*
258 CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
259 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
260 $ work, 1 )
261 CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
262 $ a( kl+i, i+1 ), lda )
263 a( kl+i, i ) = -wa
264 END IF
265*
266 IF( i.LE.min( n-1-ku, m ) ) THEN
267*
268* generate reflection to annihilate A(i,ku+i+1:n)
269*
270 wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
271 wa = sign( wn, a( i, ku+i ) )
272 IF( wn.EQ.zero ) THEN
273 tau = zero
274 ELSE
275 wb = a( i, ku+i ) + wa
276 CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
277 a( i, ku+i ) = one
278 tau = wb / wa
279 END IF
280*
281* apply reflection to A(i+1:m,ku+i:n) from the right
282*
283 CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
284 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
285 $ work, 1 )
286 CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
287 $ lda, a( i+1, ku+i ), lda )
288 a( i, ku+i ) = -wa
289 END IF
290 ELSE
291*
292* annihilate superdiagonal elements first (necessary if
293* KU = 0)
294*
295 IF( i.LE.min( n-1-ku, m ) ) THEN
296*
297* generate reflection to annihilate A(i,ku+i+1:n)
298*
299 wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
300 wa = sign( wn, a( i, ku+i ) )
301 IF( wn.EQ.zero ) THEN
302 tau = zero
303 ELSE
304 wb = a( i, ku+i ) + wa
305 CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
306 a( i, ku+i ) = one
307 tau = wb / wa
308 END IF
309*
310* apply reflection to A(i+1:m,ku+i:n) from the right
311*
312 CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
313 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
314 $ work, 1 )
315 CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
316 $ lda, a( i+1, ku+i ), lda )
317 a( i, ku+i ) = -wa
318 END IF
319*
320 IF( i.LE.min( m-1-kl, n ) ) THEN
321*
322* generate reflection to annihilate A(kl+i+1:m,i)
323*
324 wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
325 wa = sign( wn, a( kl+i, i ) )
326 IF( wn.EQ.zero ) THEN
327 tau = zero
328 ELSE
329 wb = a( kl+i, i ) + wa
330 CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
331 a( kl+i, i ) = one
332 tau = wb / wa
333 END IF
334*
335* apply reflection to A(kl+i:m,i+1:n) from the left
336*
337 CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
338 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
339 $ work, 1 )
340 CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
341 $ a( kl+i, i+1 ), lda )
342 a( kl+i, i ) = -wa
343 END IF
344 END IF
345*
346 IF (i .LE. n) THEN
347 DO 50 j = kl + i + 1, m
348 a( j, i ) = zero
349 50 CONTINUE
350 END IF
351*
352 IF (i .LE. m) THEN
353 DO 60 j = ku + i + 1, n
354 a( i, j ) = zero
355 60 CONTINUE
356 END IF
357 70 CONTINUE
358 RETURN
359*
360* End of DLAGGE
361*
362 END
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:97
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:130
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
subroutine dlagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
DLAGGE
Definition: dlagge.f:113