SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlagge.f
Go to the documentation of this file.
1 SUBROUTINE zlagge( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO )
2*
3* -- LAPACK auxiliary test routine (version 3.1) --
4* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5* November 2006
6*
7* .. Scalar Arguments ..
8 INTEGER INFO, KL, KU, LDA, M, N
9* ..
10* .. Array Arguments ..
11 INTEGER ISEED( 4 )
12 DOUBLE PRECISION D( * )
13 COMPLEX*16 A( LDA, * ), WORK( * )
14* ..
15*
16* Purpose
17* =======
18*
19* ZLAGGE generates a complex general m by n matrix A, by pre- and post-
20* multiplying a real diagonal matrix D with random unitary matrices:
21* A = U*D*V. The lower and upper bandwidths may then be reduced to
22* kl and ku by additional unitary transformations.
23*
24* Arguments
25* =========
26*
27* M (input) INTEGER
28* The number of rows of the matrix A. M >= 0.
29*
30* N (input) INTEGER
31* The number of columns of the matrix A. N >= 0.
32*
33* KL (input) INTEGER
34* The number of nonzero subdiagonals within the band of A.
35* 0 <= KL <= M-1.
36*
37* KU (input) INTEGER
38* The number of nonzero superdiagonals within the band of A.
39* 0 <= KU <= N-1.
40*
41* D (input) DOUBLE PRECISION array, dimension (min(M,N))
42* The diagonal elements of the diagonal matrix D.
43*
44* A (output) COMPLEX*16 array, dimension (LDA,N)
45* The generated m by n matrix A.
46*
47* LDA (input) INTEGER
48* The leading dimension of the array A. LDA >= M.
49*
50* ISEED (input/output) INTEGER array, dimension (4)
51* On entry, the seed of the random number generator; the array
52* elements must be between 0 and 4095, and ISEED(4) must be
53* odd.
54* On exit, the seed is updated.
55*
56* WORK (workspace) COMPLEX*16 array, dimension (M+N)
57*
58* INFO (output) INTEGER
59* = 0: successful exit
60* < 0: if INFO = -i, the i-th argument had an illegal value
61*
62* =====================================================================
63*
64* .. Parameters ..
65 COMPLEX*16 ZERO, ONE
66 parameter( zero = ( 0.0d+0, 0.0d+0 ),
67 $ one = ( 1.0d+0, 0.0d+0 ) )
68* ..
69* .. Local Scalars ..
70 INTEGER I, J
71 DOUBLE PRECISION WN
72 COMPLEX*16 TAU, WA, WB
73* ..
74* .. External Subroutines ..
75 EXTERNAL xerbla, zgemv, zgerc, zlacgv, zlarnv, zscal
76* ..
77* .. Intrinsic Functions ..
78 INTRINSIC abs, dble, max, min
79* ..
80* .. External Functions ..
81 DOUBLE PRECISION DZNRM2
82 EXTERNAL dznrm2
83* ..
84* .. Executable Statements ..
85*
86* Test the input arguments
87*
88 info = 0
89 IF( m.LT.0 ) THEN
90 info = -1
91 ELSE IF( n.LT.0 ) THEN
92 info = -2
93 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
94 info = -3
95 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
96 info = -4
97 ELSE IF( lda.LT.max( 1, m ) ) THEN
98 info = -7
99 END IF
100 IF( info.LT.0 ) THEN
101 CALL xerbla( 'ZLAGGE', -info )
102 RETURN
103 END IF
104*
105* initialize A to diagonal matrix
106*
107 DO 20 j = 1, n
108 DO 10 i = 1, m
109 a( i, j ) = zero
110 10 CONTINUE
111 20 CONTINUE
112 DO 30 i = 1, min( m, n )
113 a( i, i ) = d( i )
114 30 CONTINUE
115*
116* pre- and post-multiply A by random unitary matrices
117*
118 DO 40 i = min( m, n ), 1, -1
119 IF( i.LT.m ) THEN
120*
121* generate random reflection
122*
123 CALL zlarnv( 3, iseed, m-i+1, work )
124 wn = dznrm2( m-i+1, work, 1 )
125 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
126 IF( wn.EQ.zero ) THEN
127 tau = zero
128 ELSE
129 wb = work( 1 ) + wa
130 CALL zscal( m-i, one / wb, work( 2 ), 1 )
131 work( 1 ) = one
132 tau = dble( wb / wa )
133 END IF
134*
135* multiply A(i:m,i:n) by random reflection from the left
136*
137 CALL zgemv( 'Conjugate transpose', m-i+1, n-i+1, one,
138 $ a( i, i ), lda, work, 1, zero, work( m+1 ), 1 )
139 CALL zgerc( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
140 $ a( i, i ), lda )
141 END IF
142 IF( i.LT.n ) THEN
143*
144* generate random reflection
145*
146 CALL zlarnv( 3, iseed, n-i+1, work )
147 wn = dznrm2( n-i+1, work, 1 )
148 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
149 IF( wn.EQ.zero ) THEN
150 tau = zero
151 ELSE
152 wb = work( 1 ) + wa
153 CALL zscal( n-i, one / wb, work( 2 ), 1 )
154 work( 1 ) = one
155 tau = dble( wb / wa )
156 END IF
157*
158* multiply A(i:m,i:n) by random reflection from the right
159*
160 CALL zgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
161 $ lda, work, 1, zero, work( n+1 ), 1 )
162 CALL zgerc( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
163 $ a( i, i ), lda )
164 END IF
165 40 CONTINUE
166*
167* Reduce number of subdiagonals to KL and number of superdiagonals
168* to KU
169*
170 DO 70 i = 1, max( m-1-kl, n-1-ku )
171 IF( kl.LE.ku ) THEN
172*
173* annihilate subdiagonal elements first (necessary if KL = 0)
174*
175 IF( i.LE.min( m-1-kl, n ) ) THEN
176*
177* generate reflection to annihilate A(kl+i+1:m,i)
178*
179 wn = dznrm2( m-kl-i+1, a( kl+i, i ), 1 )
180 wa = ( wn / abs( a( kl+i, i ) ) )*a( kl+i, i )
181 IF( wn.EQ.zero ) THEN
182 tau = zero
183 ELSE
184 wb = a( kl+i, i ) + wa
185 CALL zscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
186 a( kl+i, i ) = one
187 tau = dble( wb / wa )
188 END IF
189*
190* apply reflection to A(kl+i:m,i+1:n) from the left
191*
192 CALL zgemv( 'Conjugate transpose', m-kl-i+1, n-i, one,
193 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
194 $ work, 1 )
195 CALL zgerc( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work,
196 $ 1, a( kl+i, i+1 ), lda )
197 a( kl+i, i ) = -wa
198 END IF
199*
200 IF( i.LE.min( n-1-ku, m ) ) THEN
201*
202* generate reflection to annihilate A(i,ku+i+1:n)
203*
204 wn = dznrm2( n-ku-i+1, a( i, ku+i ), lda )
205 wa = ( wn / abs( a( i, ku+i ) ) )*a( i, ku+i )
206 IF( wn.EQ.zero ) THEN
207 tau = zero
208 ELSE
209 wb = a( i, ku+i ) + wa
210 CALL zscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
211 a( i, ku+i ) = one
212 tau = dble( wb / wa )
213 END IF
214*
215* apply reflection to A(i+1:m,ku+i:n) from the right
216*
217 CALL zlacgv( n-ku-i+1, a( i, ku+i ), lda )
218 CALL zgemv( 'No transpose', m-i, n-ku-i+1, one,
219 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
220 $ work, 1 )
221 CALL zgerc( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
222 $ lda, a( i+1, ku+i ), lda )
223 a( i, ku+i ) = -wa
224 END IF
225 ELSE
226*
227* annihilate superdiagonal elements first (necessary if
228* KU = 0)
229*
230 IF( i.LE.min( n-1-ku, m ) ) THEN
231*
232* generate reflection to annihilate A(i,ku+i+1:n)
233*
234 wn = dznrm2( n-ku-i+1, a( i, ku+i ), lda )
235 wa = ( wn / abs( a( i, ku+i ) ) )*a( i, ku+i )
236 IF( wn.EQ.zero ) THEN
237 tau = zero
238 ELSE
239 wb = a( i, ku+i ) + wa
240 CALL zscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
241 a( i, ku+i ) = one
242 tau = dble( wb / wa )
243 END IF
244*
245* apply reflection to A(i+1:m,ku+i:n) from the right
246*
247 CALL zlacgv( n-ku-i+1, a( i, ku+i ), lda )
248 CALL zgemv( 'No transpose', m-i, n-ku-i+1, one,
249 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
250 $ work, 1 )
251 CALL zgerc( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
252 $ lda, a( i+1, ku+i ), lda )
253 a( i, ku+i ) = -wa
254 END IF
255*
256 IF( i.LE.min( m-1-kl, n ) ) THEN
257*
258* generate reflection to annihilate A(kl+i+1:m,i)
259*
260 wn = dznrm2( m-kl-i+1, a( kl+i, i ), 1 )
261 wa = ( wn / abs( a( kl+i, i ) ) )*a( kl+i, i )
262 IF( wn.EQ.zero ) THEN
263 tau = zero
264 ELSE
265 wb = a( kl+i, i ) + wa
266 CALL zscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
267 a( kl+i, i ) = one
268 tau = dble( wb / wa )
269 END IF
270*
271* apply reflection to A(kl+i:m,i+1:n) from the left
272*
273 CALL zgemv( 'Conjugate transpose', m-kl-i+1, n-i, one,
274 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
275 $ work, 1 )
276 CALL zgerc( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work,
277 $ 1, a( kl+i, i+1 ), lda )
278 a( kl+i, i ) = -wa
279 END IF
280 END IF
281*
282 DO 50 j = kl + i + 1, m
283 a( j, i ) = zero
284 50 CONTINUE
285*
286 DO 60 j = ku + i + 1, n
287 a( i, j ) = zero
288 60 CONTINUE
289 70 CONTINUE
290 RETURN
291*
292* End of ZLAGGE
293*
294 END
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
zlagge
subroutine zlagge(m, n, kl, ku, d, a, lda, iseed, work, info)
Definition zlagge.f:2
zlarnv
subroutine zlarnv(idist, iseed, n, x)
Definition zlarnv.f:2