LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.

Definition at line 112 of file dlagge.f.

113*
114* -- LAPACK auxiliary routine --
115* -- LAPACK is a software package provided by Univ. of Tennessee, --
116* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
117*
118* .. Scalar Arguments ..
119 INTEGER INFO, KL, KU, LDA, M, N
120* ..
121* .. Array Arguments ..
122 INTEGER ISEED( 4 )
123 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
124* ..
125*
126* =====================================================================
127*
128* .. Parameters ..
129 DOUBLE PRECISION ZERO, ONE
130 parameter( zero = 0.0d+0, one = 1.0d+0 )
131* ..
132* .. Local Scalars ..
133 INTEGER I, J
134 DOUBLE PRECISION TAU, WA, WB, WN
135* ..
136* .. External Subroutines ..
137 EXTERNAL dgemv, dger, dlarnv, dscal, xerbla
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC max, min, sign
141* ..
142* .. External Functions ..
143 DOUBLE PRECISION DNRM2
144 EXTERNAL dnrm2
145* ..
146* .. Executable Statements ..
147*
148* Test the input arguments
149*
150 info = 0
151 IF( m.LT.0 ) THEN
152 info = -1
153 ELSE IF( n.LT.0 ) THEN
154 info = -2
155 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
156 info = -3
157 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
158 info = -4
159 ELSE IF( lda.LT.max( 1, m ) ) THEN
160 info = -7
161 END IF
162 IF( info.LT.0 ) THEN
163 CALL xerbla( 'DLAGGE', -info )
164 RETURN
165 END IF
166*
167* initialize A to diagonal matrix
168*
169 DO 20 j = 1, n
170 DO 10 i = 1, m
171 a( i, j ) = zero
172 10 CONTINUE
173 20 CONTINUE
174 DO 30 i = 1, min( m, n )
175 a( i, i ) = d( i )
176 30 CONTINUE
177*
178* Quick exit if the user wants a diagonal matrix
179*
180 IF(( kl .EQ. 0 ).AND.( ku .EQ. 0)) RETURN
181*
182* pre- and post-multiply A by random orthogonal matrices
183*
184 DO 40 i = min( m, n ), 1, -1
185 IF( i.LT.m ) THEN
186*
187* generate random reflection
188*
189 CALL dlarnv( 3, iseed, m-i+1, work )
190 wn = dnrm2( m-i+1, work, 1 )
191 wa = sign( wn, work( 1 ) )
192 IF( wn.EQ.zero ) THEN
193 tau = zero
194 ELSE
195 wb = work( 1 ) + wa
196 CALL dscal( m-i, one / wb, work( 2 ), 1 )
197 work( 1 ) = one
198 tau = wb / wa
199 END IF
200*
201* multiply A(i:m,i:n) by random reflection from the left
202*
203 CALL dgemv( 'Transpose', m-i+1, n-i+1, one, a( i, i ), lda,
204 $ work, 1, zero, work( m+1 ), 1 )
205 CALL dger( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
206 $ a( i, i ), lda )
207 END IF
208 IF( i.LT.n ) THEN
209*
210* generate random reflection
211*
212 CALL dlarnv( 3, iseed, n-i+1, work )
213 wn = dnrm2( n-i+1, work, 1 )
214 wa = sign( wn, work( 1 ) )
215 IF( wn.EQ.zero ) THEN
216 tau = zero
217 ELSE
218 wb = work( 1 ) + wa
219 CALL dscal( n-i, one / wb, work( 2 ), 1 )
220 work( 1 ) = one
221 tau = wb / wa
222 END IF
223*
224* multiply A(i:m,i:n) by random reflection from the right
225*
226 CALL dgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
227 $ lda, work, 1, zero, work( n+1 ), 1 )
228 CALL dger( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
229 $ a( i, i ), lda )
230 END IF
231 40 CONTINUE
232*
233* Reduce number of subdiagonals to KL and number of superdiagonals
234* to KU
235*
236 DO 70 i = 1, max( m-1-kl, n-1-ku )
237 IF( kl.LE.ku ) THEN
238*
239* annihilate subdiagonal elements first (necessary if KL = 0)
240*
241 IF( i.LE.min( m-1-kl, n ) ) THEN
242*
243* generate reflection to annihilate A(kl+i+1:m,i)
244*
245 wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
246 wa = sign( wn, a( kl+i, i ) )
247 IF( wn.EQ.zero ) THEN
248 tau = zero
249 ELSE
250 wb = a( kl+i, i ) + wa
251 CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
252 a( kl+i, i ) = one
253 tau = wb / wa
254 END IF
255*
256* apply reflection to A(kl+i:m,i+1:n) from the left
257*
258 CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
259 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
260 $ work, 1 )
261 CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
262 $ a( kl+i, i+1 ), lda )
263 a( kl+i, i ) = -wa
264 END IF
265*
266 IF( i.LE.min( n-1-ku, m ) ) THEN
267*
268* generate reflection to annihilate A(i,ku+i+1:n)
269*
270 wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
271 wa = sign( wn, a( i, ku+i ) )
272 IF( wn.EQ.zero ) THEN
273 tau = zero
274 ELSE
275 wb = a( i, ku+i ) + wa
276 CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
277 a( i, ku+i ) = one
278 tau = wb / wa
279 END IF
280*
281* apply reflection to A(i+1:m,ku+i:n) from the right
282*
283 CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
284 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
285 $ work, 1 )
286 CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
287 $ lda, a( i+1, ku+i ), lda )
288 a( i, ku+i ) = -wa
289 END IF
290 ELSE
291*
292* annihilate superdiagonal elements first (necessary if
293* KU = 0)
294*
295 IF( i.LE.min( n-1-ku, m ) ) THEN
296*
297* generate reflection to annihilate A(i,ku+i+1:n)
298*
299 wn = dnrm2( n-ku-i+1, a( i, ku+i ), lda )
300 wa = sign( wn, a( i, ku+i ) )
301 IF( wn.EQ.zero ) THEN
302 tau = zero
303 ELSE
304 wb = a( i, ku+i ) + wa
305 CALL dscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
306 a( i, ku+i ) = one
307 tau = wb / wa
308 END IF
309*
310* apply reflection to A(i+1:m,ku+i:n) from the right
311*
312 CALL dgemv( 'No transpose', m-i, n-ku-i+1, one,
313 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
314 $ work, 1 )
315 CALL dger( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
316 $ lda, a( i+1, ku+i ), lda )
317 a( i, ku+i ) = -wa
318 END IF
319*
320 IF( i.LE.min( m-1-kl, n ) ) THEN
321*
322* generate reflection to annihilate A(kl+i+1:m,i)
323*
324 wn = dnrm2( m-kl-i+1, a( kl+i, i ), 1 )
325 wa = sign( wn, a( kl+i, i ) )
326 IF( wn.EQ.zero ) THEN
327 tau = zero
328 ELSE
329 wb = a( kl+i, i ) + wa
330 CALL dscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
331 a( kl+i, i ) = one
332 tau = wb / wa
333 END IF
334*
335* apply reflection to A(kl+i:m,i+1:n) from the left
336*
337 CALL dgemv( 'Transpose', m-kl-i+1, n-i, one,
338 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
339 $ work, 1 )
340 CALL dger( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work, 1,
341 $ a( kl+i, i+1 ), lda )
342 a( kl+i, i ) = -wa
343 END IF
344 END IF
345*
346 IF (i .LE. n) THEN
347 DO 50 j = kl + i + 1, m
348 a( j, i ) = zero
349 50 CONTINUE
350 END IF
351*
352 IF (i .LE. m) THEN
353 DO 60 j = ku + i + 1, n
354 a( i, j ) = zero
355 60 CONTINUE
356 END IF
357 70 CONTINUE
358 RETURN
359*
360* End of DLAGGE
361*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
dlarnv
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:95
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: