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