LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERC
Definition: cgerc.f:130
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: clarnv.f:99
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
Here is the call graph for this function:
Here is the caller graph for this function: