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