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
369