LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlagge.f
Go to the documentation of this file.
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