ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
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
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
min
#define min(A, B)
Definition: pcgemr.c:181