LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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)
XERBLA
Definition: xerbla.f:60
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERC
Definition: zgerc.f:130
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zlagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
ZLAGGE
Definition: zlagge.f:114
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: zlarnv.f:99
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74