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

◆ clagge()

subroutine clagge ( integer  m,
integer  n,
integer  kl,
integer  ku,
real, dimension( * )  d,
complex, dimension( lda, * )  a,
integer  lda,
integer, dimension( 4 )  iseed,
complex, dimension( * )  work,
integer  info 
)

CLAGGE

Purpose:
 CLAGGE generates a complex general m by n matrix A, by pre- and post-
 multiplying a real diagonal matrix D with random unitary matrices:
 A = U*D*V. The lower and upper bandwidths may then be reduced to
 kl and ku by additional unitary 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 REAL array, dimension (min(M,N))
          The diagonal elements of the diagonal matrix D.
[out]A
          A is COMPLEX 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 COMPLEX 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 113 of file clagge.f.

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 REAL D( * )
125 COMPLEX A( LDA, * ), WORK( * )
126* ..
127*
128* =====================================================================
129*
130* .. Parameters ..
131 COMPLEX ZERO, ONE
132 parameter( zero = ( 0.0e+0, 0.0e+0 ),
133 $ one = ( 1.0e+0, 0.0e+0 ) )
134* ..
135* .. Local Scalars ..
136 INTEGER I, J
137 REAL WN
138 COMPLEX TAU, WA, WB
139* ..
140* .. External Subroutines ..
141 EXTERNAL cgemv, cgerc, clacgv, clarnv, cscal, xerbla
142* ..
143* .. Intrinsic Functions ..
144 INTRINSIC abs, max, min, real
145* ..
146* .. External Functions ..
147 REAL SCNRM2
148 EXTERNAL scnrm2
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( 'CLAGGE', -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 clarnv( 3, iseed, m-i+1, work )
194 wn = scnrm2( 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 cscal( m-i, one / wb, work( 2 ), 1 )
201 work( 1 ) = one
202 tau = real( wb / wa )
203 END IF
204*
205* multiply A(i:m,i:n) by random reflection from the left
206*
207 CALL cgemv( 'Conjugate transpose', m-i+1, n-i+1, one,
208 $ a( i, i ), lda, work, 1, zero, work( m+1 ), 1 )
209 CALL cgerc( 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 clarnv( 3, iseed, n-i+1, work )
217 wn = scnrm2( 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 cscal( n-i, one / wb, work( 2 ), 1 )
224 work( 1 ) = one
225 tau = real( wb / wa )
226 END IF
227*
228* multiply A(i:m,i:n) by random reflection from the right
229*
230 CALL cgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
231 $ lda, work, 1, zero, work( n+1 ), 1 )
232 CALL cgerc( 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 = scnrm2( 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 cscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
256 a( kl+i, i ) = one
257 tau = real( wb / wa )
258 END IF
259*
260* apply reflection to A(kl+i:m,i+1:n) from the left
261*
262 CALL cgemv( '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 cgerc( 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 = scnrm2( 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 cscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
281 a( i, ku+i ) = one
282 tau = real( wb / wa )
283 END IF
284*
285* apply reflection to A(i+1:m,ku+i:n) from the right
286*
287 CALL clacgv( n-ku-i+1, a( i, ku+i ), lda )
288 CALL cgemv( '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 cgerc( 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 = scnrm2( 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 cscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
311 a( i, ku+i ) = one
312 tau = real( wb / wa )
313 END IF
314*
315* apply reflection to A(i+1:m,ku+i:n) from the right
316*
317 CALL clacgv( n-ku-i+1, a( i, ku+i ), lda )
318 CALL cgemv( '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 cgerc( 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 = scnrm2( 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 cscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
337 a( kl+i, i ) = one
338 tau = real( wb / wa )
339 END IF
340*
341* apply reflection to A(kl+i:m,i+1:n) from the left
342*
343 CALL cgemv( '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 cgerc( 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 CLAGGE
367*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
cgerc
subroutine cgerc(m, n, alpha, x, incx, y, incy, a, lda)
CGERC
Definition cgerc.f:130
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
clarnv
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:99
scnrm2
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: