LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dlagge()

subroutine dlagge ( integer  M,
integer  N,
integer  KL,
integer  KU,
double precision, dimension( * )  D,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( 4 )  ISEED,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLAGGE

Purpose:
 DLAGGE generates a real general m by n matrix A, by pre- and post-
 multiplying a real diagonal matrix D with random orthogonal matrices:
 A = U*D*V. The lower and upper bandwidths may then be reduced to
 kl and ku by additional orthogonal transformations.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]KL
          KL is INTEGER
          The number of nonzero subdiagonals within the band of A.
          0 <= KL <= M-1.
[in]KU
          KU is INTEGER
          The number of nonzero superdiagonals within the band of A.
          0 <= KU <= N-1.
[in]D
          D is DOUBLE PRECISION array, dimension (min(M,N))
          The diagonal elements of the diagonal matrix D.
[out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The generated m by n matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= M.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry, the seed of the random number generator; the array
          elements must be between 0 and 4095, and ISEED(4) must be
          odd.
          On exit, the seed is updated.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (M+N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 115 of file dlagge.f.

115 *
116 * -- LAPACK auxiliary routine (version 3.7.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 * December 2016
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 * Quick exit if the user wants a diagonal matrix
182 *
183  IF(( kl .EQ. 0 ).AND.( ku .EQ. 0)) RETURN
184 *
185 * pre- and post-multiply A by random orthogonal matrices
186 *
187  DO 40 i = min( m, n ), 1, -1
188  IF( i.LT.m ) THEN
189 *
190 * generate random reflection
191 *
192  CALL dlarnv( 3, iseed, m-i+1, work )
193  wn = dnrm2( m-i+1, work, 1 )
194  wa = sign( wn, work( 1 ) )
195  IF( wn.EQ.zero ) THEN
196  tau = zero
197  ELSE
198  wb = work( 1 ) + wa
199  CALL dscal( m-i, one / wb, work( 2 ), 1 )
200  work( 1 ) = one
201  tau = wb / wa
202  END IF
203 *
204 * multiply A(i:m,i:n) by random reflection from the left
205 *
206  CALL dgemv( 'Transpose', m-i+1, n-i+1, one, a( i, i ), lda,
207  $ work, 1, zero, work( m+1 ), 1 )
208  CALL dger( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
209  $ a( i, i ), lda )
210  END IF
211  IF( i.LT.n ) THEN
212 *
213 * generate random reflection
214 *
215  CALL dlarnv( 3, iseed, n-i+1, work )
216  wn = dnrm2( n-i+1, work, 1 )
217  wa = sign( wn, work( 1 ) )
218  IF( wn.EQ.zero ) THEN
219  tau = zero
220  ELSE
221  wb = work( 1 ) + wa
222  CALL dscal( n-i, one / wb, work( 2 ), 1 )
223  work( 1 ) = one
224  tau = wb / wa
225  END IF
226 *
227 * multiply A(i:m,i:n) by random reflection from the right
228 *
229  CALL dgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
230  $ lda, work, 1, zero, work( n+1 ), 1 )
231  CALL dger( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
232  $ a( i, i ), lda )
233  END IF
234  40 CONTINUE
235 *
236 * Reduce number of subdiagonals to KL and number of superdiagonals
237 * to KU
238 *
239  DO 70 i = 1, max( m-1-kl, n-1-ku )
240  IF( kl.LE.ku ) THEN
241 *
242 * annihilate subdiagonal elements first (necessary if KL = 0)
243 *
244  IF( i.LE.min( m-1-kl, n ) ) THEN
245 *
246 * generate reflection to annihilate A(kl+i+1:m,i)
247 *
248  wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
249  wa = sign( wn, a( kl+i, i ) )
250  IF( wn.EQ.zero ) THEN
251  tau = zero
252  ELSE
253  wb = a( kl+i, i ) + wa
254  CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
255  a( kl+i, i ) = one
256  tau = wb / wa
257  END IF
258 *
259 * apply reflection to A(kl+i:m,i+1:n) from the left
260 *
261  CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
262  $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
263  $ work, 1 )
264  CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
265  $ a( kl+i, i+1 ), lda )
266  a( kl+i, i ) = -wa
267  END IF
268 *
269  IF( i.LE.min( n-1-ku, m ) ) THEN
270 *
271 * generate reflection to annihilate A(i,ku+i+1:n)
272 *
273  wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
274  wa = sign( wn, a( i, ku+i ) )
275  IF( wn.EQ.zero ) THEN
276  tau = zero
277  ELSE
278  wb = a( i, ku+i ) + wa
279  CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
280  a( i, ku+i ) = one
281  tau = wb / wa
282  END IF
283 *
284 * apply reflection to A(i+1:m,ku+i:n) from the right
285 *
286  CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
287  $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
288  $ work, 1 )
289  CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
290  $ lda, a( i+1, ku+i ), lda )
291  a( i, ku+i ) = -wa
292  END IF
293  ELSE
294 *
295 * annihilate superdiagonal elements first (necessary if
296 * KU = 0)
297 *
298  IF( i.LE.min( n-1-ku, m ) ) THEN
299 *
300 * generate reflection to annihilate A(i,ku+i+1:n)
301 *
302  wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
303  wa = sign( wn, a( i, ku+i ) )
304  IF( wn.EQ.zero ) THEN
305  tau = zero
306  ELSE
307  wb = a( i, ku+i ) + wa
308  CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
309  a( i, ku+i ) = one
310  tau = wb / wa
311  END IF
312 *
313 * apply reflection to A(i+1:m,ku+i:n) from the right
314 *
315  CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
316  $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
317  $ work, 1 )
318  CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
319  $ lda, a( i+1, ku+i ), lda )
320  a( i, ku+i ) = -wa
321  END IF
322 *
323  IF( i.LE.min( m-1-kl, n ) ) THEN
324 *
325 * generate reflection to annihilate A(kl+i+1:m,i)
326 *
327  wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
328  wa = sign( wn, a( kl+i, i ) )
329  IF( wn.EQ.zero ) THEN
330  tau = zero
331  ELSE
332  wb = a( kl+i, i ) + wa
333  CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
334  a( kl+i, i ) = one
335  tau = wb / wa
336  END IF
337 *
338 * apply reflection to A(kl+i:m,i+1:n) from the left
339 *
340  CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
341  $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
342  $ work, 1 )
343  CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
344  $ a( kl+i, i+1 ), lda )
345  a( kl+i, i ) = -wa
346  END IF
347  END IF
348 *
349  IF (i .LE. n) THEN
350  DO 50 j = kl + i + 1, m
351  a( j, i ) = zero
352  50 CONTINUE
353  END IF
354 *
355  IF (i .LE. m) THEN
356  DO 60 j = ku + i + 1, n
357  a( i, j ) = zero
358  60 CONTINUE
359  END IF
360  70 CONTINUE
361  RETURN
362 *
363 * End of DLAGGE
364 *
Here is the call graph for this function:
Here is the caller graph for this function:
dger
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
dnrm2
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
dlarnv
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81