LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
clatms.f
Go to the documentation of this file.
1 *> \brief \b CLATMS
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 CLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
12 * KL, KU, PACK, A, LDA, WORK, INFO )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER DIST, PACK, SYM
16 * INTEGER INFO, KL, KU, LDA, M, MODE, N
17 * REAL COND, DMAX
18 * ..
19 * .. Array Arguments ..
20 * INTEGER ISEED( 4 )
21 * REAL D( * )
22 * COMPLEX A( LDA, * ), WORK( * )
23 * ..
24 *
25 *
26 *> \par Purpose:
27 * =============
28 *>
29 *> \verbatim
30 *>
31 *> CLATMS generates random matrices with specified singular values
32 *> (or hermitian with specified eigenvalues)
33 *> for testing LAPACK programs.
34 *>
35 *> CLATMS operates by applying the following sequence of
36 *> operations:
37 *>
38 *> Set the diagonal to D, where D may be input or
39 *> computed according to MODE, COND, DMAX, and SYM
40 *> as described below.
41 *>
42 *> Generate a matrix with the appropriate band structure, by one
43 *> of two methods:
44 *>
45 *> Method A:
46 *> Generate a dense M x N matrix by multiplying D on the left
47 *> and the right by random unitary matrices, then:
48 *>
49 *> Reduce the bandwidth according to KL and KU, using
50 *> Householder transformations.
51 *>
52 *> Method B:
53 *> Convert the bandwidth-0 (i.e., diagonal) matrix to a
54 *> bandwidth-1 matrix using Givens rotations, "chasing"
55 *> out-of-band elements back, much as in QR; then convert
56 *> the bandwidth-1 to a bandwidth-2 matrix, etc. Note
57 *> that for reasonably small bandwidths (relative to M and
58 *> N) this requires less storage, as a dense matrix is not
59 *> generated. Also, for hermitian or symmetric matrices,
60 *> only one triangle is generated.
61 *>
62 *> Method A is chosen if the bandwidth is a large fraction of the
63 *> order of the matrix, and LDA is at least M (so a dense
64 *> matrix can be stored.) Method B is chosen if the bandwidth
65 *> is small (< 1/2 N for hermitian or symmetric, < .3 N+M for
66 *> non-symmetric), or LDA is less than M and not less than the
67 *> bandwidth.
68 *>
69 *> Pack the matrix if desired. Options specified by PACK are:
70 *> no packing
71 *> zero out upper half (if hermitian)
72 *> zero out lower half (if hermitian)
73 *> store the upper half columnwise (if hermitian or upper
74 *> triangular)
75 *> store the lower half columnwise (if hermitian or lower
76 *> triangular)
77 *> store the lower triangle in banded format (if hermitian or
78 *> lower triangular)
79 *> store the upper triangle in banded format (if hermitian or
80 *> upper triangular)
81 *> store the entire matrix in banded format
82 *> If Method B is chosen, and band format is specified, then the
83 *> matrix will be generated in the band format, so no repacking
84 *> will be necessary.
85 *> \endverbatim
86 *
87 * Arguments:
88 * ==========
89 *
90 *> \param[in] M
91 *> \verbatim
92 *> M is INTEGER
93 *> The number of rows of A. Not modified.
94 *> \endverbatim
95 *>
96 *> \param[in] N
97 *> \verbatim
98 *> N is INTEGER
99 *> The number of columns of A. N must equal M if the matrix
100 *> is symmetric or hermitian (i.e., if SYM is not 'N')
101 *> Not modified.
102 *> \endverbatim
103 *>
104 *> \param[in] DIST
105 *> \verbatim
106 *> DIST is CHARACTER*1
107 *> On entry, DIST specifies the type of distribution to be used
108 *> to generate the random eigen-/singular values.
109 *> 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
110 *> 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
111 *> 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
112 *> Not modified.
113 *> \endverbatim
114 *>
115 *> \param[in,out] ISEED
116 *> \verbatim
117 *> ISEED is INTEGER array, dimension ( 4 )
118 *> On entry ISEED specifies the seed of the random number
119 *> generator. They should lie between 0 and 4095 inclusive,
120 *> and ISEED(4) should be odd. The random number generator
121 *> uses a linear congruential sequence limited to small
122 *> integers, and so should produce machine independent
123 *> random numbers. The values of ISEED are changed on
124 *> exit, and can be used in the next call to CLATMS
125 *> to continue the same random number sequence.
126 *> Changed on exit.
127 *> \endverbatim
128 *>
129 *> \param[in] SYM
130 *> \verbatim
131 *> SYM is CHARACTER*1
132 *> If SYM='H', the generated matrix is hermitian, with
133 *> eigenvalues specified by D, COND, MODE, and DMAX; they
134 *> may be positive, negative, or zero.
135 *> If SYM='P', the generated matrix is hermitian, with
136 *> eigenvalues (= singular values) specified by D, COND,
137 *> MODE, and DMAX; they will not be negative.
138 *> If SYM='N', the generated matrix is nonsymmetric, with
139 *> singular values specified by D, COND, MODE, and DMAX;
140 *> they will not be negative.
141 *> If SYM='S', the generated matrix is (complex) symmetric,
142 *> with singular values specified by D, COND, MODE, and
143 *> DMAX; they will not be negative.
144 *> Not modified.
145 *> \endverbatim
146 *>
147 *> \param[in,out] D
148 *> \verbatim
149 *> D is REAL array, dimension ( MIN( M, N ) )
150 *> This array is used to specify the singular values or
151 *> eigenvalues of A (see SYM, above.) If MODE=0, then D is
152 *> assumed to contain the singular/eigenvalues, otherwise
153 *> they will be computed according to MODE, COND, and DMAX,
154 *> and placed in D.
155 *> Modified if MODE is nonzero.
156 *> \endverbatim
157 *>
158 *> \param[in] MODE
159 *> \verbatim
160 *> MODE is INTEGER
161 *> On entry this describes how the singular/eigenvalues are to
162 *> be specified:
163 *> MODE = 0 means use D as input
164 *> MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
165 *> MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
166 *> MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
167 *> MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
168 *> MODE = 5 sets D to random numbers in the range
169 *> ( 1/COND , 1 ) such that their logarithms
170 *> are uniformly distributed.
171 *> MODE = 6 set D to random numbers from same distribution
172 *> as the rest of the matrix.
173 *> MODE < 0 has the same meaning as ABS(MODE), except that
174 *> the order of the elements of D is reversed.
175 *> Thus if MODE is positive, D has entries ranging from
176 *> 1 to 1/COND, if negative, from 1/COND to 1,
177 *> If SYM='H', and MODE is neither 0, 6, nor -6, then
178 *> the elements of D will also be multiplied by a random
179 *> sign (i.e., +1 or -1.)
180 *> Not modified.
181 *> \endverbatim
182 *>
183 *> \param[in] COND
184 *> \verbatim
185 *> COND is REAL
186 *> On entry, this is used as described under MODE above.
187 *> If used, it must be >= 1. Not modified.
188 *> \endverbatim
189 *>
190 *> \param[in] DMAX
191 *> \verbatim
192 *> DMAX is REAL
193 *> If MODE is neither -6, 0 nor 6, the contents of D, as
194 *> computed according to MODE and COND, will be scaled by
195 *> DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
196 *> singular value (which is to say the norm) will be abs(DMAX).
197 *> Note that DMAX need not be positive: if DMAX is negative
198 *> (or zero), D will be scaled by a negative number (or zero).
199 *> Not modified.
200 *> \endverbatim
201 *>
202 *> \param[in] KL
203 *> \verbatim
204 *> KL is INTEGER
205 *> This specifies the lower bandwidth of the matrix. For
206 *> example, KL=0 implies upper triangular, KL=1 implies upper
207 *> Hessenberg, and KL being at least M-1 means that the matrix
208 *> has full lower bandwidth. KL must equal KU if the matrix
209 *> is symmetric or hermitian.
210 *> Not modified.
211 *> \endverbatim
212 *>
213 *> \param[in] KU
214 *> \verbatim
215 *> KU is INTEGER
216 *> This specifies the upper bandwidth of the matrix. For
217 *> example, KU=0 implies lower triangular, KU=1 implies lower
218 *> Hessenberg, and KU being at least N-1 means that the matrix
219 *> has full upper bandwidth. KL must equal KU if the matrix
220 *> is symmetric or hermitian.
221 *> Not modified.
222 *> \endverbatim
223 *>
224 *> \param[in] PACK
225 *> \verbatim
226 *> PACK is CHARACTER*1
227 *> This specifies packing of matrix as follows:
228 *> 'N' => no packing
229 *> 'U' => zero out all subdiagonal entries (if symmetric
230 *> or hermitian)
231 *> 'L' => zero out all superdiagonal entries (if symmetric
232 *> or hermitian)
233 *> 'C' => store the upper triangle columnwise (only if the
234 *> matrix is symmetric, hermitian, or upper triangular)
235 *> 'R' => store the lower triangle columnwise (only if the
236 *> matrix is symmetric, hermitian, or lower triangular)
237 *> 'B' => store the lower triangle in band storage scheme
238 *> (only if the matrix is symmetric, hermitian, or
239 *> lower triangular)
240 *> 'Q' => store the upper triangle in band storage scheme
241 *> (only if the matrix is symmetric, hermitian, or
242 *> upper triangular)
243 *> 'Z' => store the entire matrix in band storage scheme
244 *> (pivoting can be provided for by using this
245 *> option to store A in the trailing rows of
246 *> the allocated storage)
247 *>
248 *> Using these options, the various LAPACK packed and banded
249 *> storage schemes can be obtained:
250 *> GB - use 'Z'
251 *> PB, SB, HB, or TB - use 'B' or 'Q'
252 *> PP, SP, HB, or TP - use 'C' or 'R'
253 *>
254 *> If two calls to CLATMS differ only in the PACK parameter,
255 *> they will generate mathematically equivalent matrices.
256 *> Not modified.
257 *> \endverbatim
258 *>
259 *> \param[in,out] A
260 *> \verbatim
261 *> A is COMPLEX array, dimension ( LDA, N )
262 *> On exit A is the desired test matrix. A is first generated
263 *> in full (unpacked) form, and then packed, if so specified
264 *> by PACK. Thus, the first M elements of the first N
265 *> columns will always be modified. If PACK specifies a
266 *> packed or banded storage scheme, all LDA elements of the
267 *> first N columns will be modified; the elements of the
268 *> array which do not correspond to elements of the generated
269 *> matrix are set to zero.
270 *> Modified.
271 *> \endverbatim
272 *>
273 *> \param[in] LDA
274 *> \verbatim
275 *> LDA is INTEGER
276 *> LDA specifies the first dimension of A as declared in the
277 *> calling program. If PACK='N', 'U', 'L', 'C', or 'R', then
278 *> LDA must be at least M. If PACK='B' or 'Q', then LDA must
279 *> be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
280 *> If PACK='Z', LDA must be large enough to hold the packed
281 *> array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
282 *> Not modified.
283 *> \endverbatim
284 *>
285 *> \param[out] WORK
286 *> \verbatim
287 *> WORK is COMPLEX array, dimension ( 3*MAX( N, M ) )
288 *> Workspace.
289 *> Modified.
290 *> \endverbatim
291 *>
292 *> \param[out] INFO
293 *> \verbatim
294 *> INFO is INTEGER
295 *> Error code. On exit, INFO will be set to one of the
296 *> following values:
297 *> 0 => normal return
298 *> -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
299 *> -2 => N negative
300 *> -3 => DIST illegal string
301 *> -5 => SYM illegal string
302 *> -7 => MODE not in range -6 to 6
303 *> -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
304 *> -10 => KL negative
305 *> -11 => KU negative, or SYM is not 'N' and KU is not equal to
306 *> KL
307 *> -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
308 *> or PACK='C' or 'Q' and SYM='N' and KL is not zero;
309 *> or PACK='R' or 'B' and SYM='N' and KU is not zero;
310 *> or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
311 *> N.
312 *> -14 => LDA is less than M, or PACK='Z' and LDA is less than
313 *> MIN(KU,N-1) + MIN(KL,M-1) + 1.
314 *> 1 => Error return from SLATM1
315 *> 2 => Cannot scale to DMAX (max. sing. value is 0)
316 *> 3 => Error return from CLAGGE, CLAGHE or CLAGSY
317 *> \endverbatim
318 *
319 * Authors:
320 * ========
321 *
322 *> \author Univ. of Tennessee
323 *> \author Univ. of California Berkeley
324 *> \author Univ. of Colorado Denver
325 *> \author NAG Ltd.
326 *
327 *> \ingroup complex_matgen
328 *
329 * =====================================================================
330  SUBROUTINE clatms( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
331  $ KL, KU, PACK, A, LDA, WORK, INFO )
332 *
333 * -- LAPACK computational routine --
334 * -- LAPACK is a software package provided by Univ. of Tennessee, --
335 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
336 *
337 * .. Scalar Arguments ..
338  CHARACTER DIST, PACK, SYM
339  INTEGER INFO, KL, KU, LDA, M, MODE, N
340  REAL COND, DMAX
341 * ..
342 * .. Array Arguments ..
343  INTEGER ISEED( 4 )
344  REAL D( * )
345  COMPLEX A( LDA, * ), WORK( * )
346 * ..
347 *
348 * =====================================================================
349 *
350 * .. Parameters ..
351  REAL ZERO
352  parameter( zero = 0.0e+0 )
353  REAL ONE
354  parameter( one = 1.0e+0 )
355  COMPLEX CZERO
356  parameter( czero = ( 0.0e+0, 0.0e+0 ) )
357  REAL TWOPI
358  parameter( twopi = 6.28318530717958647692528676655900576839e+0 )
359 * ..
360 * .. Local Scalars ..
361  LOGICAL CSYM, GIVENS, ILEXTR, ILTEMP, TOPDWN
362  INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
363  $ ioffg, ioffst, ipack, ipackg, ir, ir1, ir2,
364  $ irow, irsign, iskew, isym, isympk, j, jc, jch,
365  $ jkl, jku, jr, k, llb, minlda, mnmin, mr, nc,
366  $ uub
367  REAL ALPHA, ANGLE, REALC, TEMP
368  COMPLEX C, CT, CTEMP, DUMMY, EXTRA, S, ST
369 * ..
370 * .. External Functions ..
371  LOGICAL LSAME
372  REAL SLARND
373  COMPLEX CLARND
374  EXTERNAL lsame, slarnd, clarnd
375 * ..
376 * .. External Subroutines ..
377  EXTERNAL clagge, claghe, clagsy, clarot, clartg, claset,
378  $ slatm1, sscal, xerbla
379 * ..
380 * .. Intrinsic Functions ..
381  INTRINSIC abs, cmplx, conjg, cos, max, min, mod, real,
382  $ sin
383 * ..
384 * .. Executable Statements ..
385 *
386 * 1) Decode and Test the input parameters.
387 * Initialize flags & seed.
388 *
389  info = 0
390 *
391 * Quick return if possible
392 *
393  IF( m.EQ.0 .OR. n.EQ.0 )
394  $ RETURN
395 *
396 * Decode DIST
397 *
398  IF( lsame( dist, 'U' ) ) THEN
399  idist = 1
400  ELSE IF( lsame( dist, 'S' ) ) THEN
401  idist = 2
402  ELSE IF( lsame( dist, 'N' ) ) THEN
403  idist = 3
404  ELSE
405  idist = -1
406  END IF
407 *
408 * Decode SYM
409 *
410  IF( lsame( sym, 'N' ) ) THEN
411  isym = 1
412  irsign = 0
413  csym = .false.
414  ELSE IF( lsame( sym, 'P' ) ) THEN
415  isym = 2
416  irsign = 0
417  csym = .false.
418  ELSE IF( lsame( sym, 'S' ) ) THEN
419  isym = 2
420  irsign = 0
421  csym = .true.
422  ELSE IF( lsame( sym, 'H' ) ) THEN
423  isym = 2
424  irsign = 1
425  csym = .false.
426  ELSE
427  isym = -1
428  END IF
429 *
430 * Decode PACK
431 *
432  isympk = 0
433  IF( lsame( pack, 'N' ) ) THEN
434  ipack = 0
435  ELSE IF( lsame( pack, 'U' ) ) THEN
436  ipack = 1
437  isympk = 1
438  ELSE IF( lsame( pack, 'L' ) ) THEN
439  ipack = 2
440  isympk = 1
441  ELSE IF( lsame( pack, 'C' ) ) THEN
442  ipack = 3
443  isympk = 2
444  ELSE IF( lsame( pack, 'R' ) ) THEN
445  ipack = 4
446  isympk = 3
447  ELSE IF( lsame( pack, 'B' ) ) THEN
448  ipack = 5
449  isympk = 3
450  ELSE IF( lsame( pack, 'Q' ) ) THEN
451  ipack = 6
452  isympk = 2
453  ELSE IF( lsame( pack, 'Z' ) ) THEN
454  ipack = 7
455  ELSE
456  ipack = -1
457  END IF
458 *
459 * Set certain internal parameters
460 *
461  mnmin = min( m, n )
462  llb = min( kl, m-1 )
463  uub = min( ku, n-1 )
464  mr = min( m, n+llb )
465  nc = min( n, m+uub )
466 *
467  IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
468  minlda = uub + 1
469  ELSE IF( ipack.EQ.7 ) THEN
470  minlda = llb + uub + 1
471  ELSE
472  minlda = m
473  END IF
474 *
475 * Use Givens rotation method if bandwidth small enough,
476 * or if LDA is too small to store the matrix unpacked.
477 *
478  givens = .false.
479  IF( isym.EQ.1 ) THEN
480  IF( real( llb+uub ).LT.0.3*real( max( 1, mr+nc ) ) )
481  $ givens = .true.
482  ELSE
483  IF( 2*llb.LT.m )
484  $ givens = .true.
485  END IF
486  IF( lda.LT.m .AND. lda.GE.minlda )
487  $ givens = .true.
488 *
489 * Set INFO if an error
490 *
491  IF( m.LT.0 ) THEN
492  info = -1
493  ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
494  info = -1
495  ELSE IF( n.LT.0 ) THEN
496  info = -2
497  ELSE IF( idist.EQ.-1 ) THEN
498  info = -3
499  ELSE IF( isym.EQ.-1 ) THEN
500  info = -5
501  ELSE IF( abs( mode ).GT.6 ) THEN
502  info = -7
503  ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
504  $ THEN
505  info = -8
506  ELSE IF( kl.LT.0 ) THEN
507  info = -10
508  ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
509  info = -11
510  ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
511  $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
512  $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
513  $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
514  info = -12
515  ELSE IF( lda.LT.max( 1, minlda ) ) THEN
516  info = -14
517  END IF
518 *
519  IF( info.NE.0 ) THEN
520  CALL xerbla( 'CLATMS', -info )
521  RETURN
522  END IF
523 *
524 * Initialize random number generator
525 *
526  DO 10 i = 1, 4
527  iseed( i ) = mod( abs( iseed( i ) ), 4096 )
528  10 CONTINUE
529 *
530  IF( mod( iseed( 4 ), 2 ).NE.1 )
531  $ iseed( 4 ) = iseed( 4 ) + 1
532 *
533 * 2) Set up D if indicated.
534 *
535 * Compute D according to COND and MODE
536 *
537  CALL slatm1( mode, cond, irsign, idist, iseed, d, mnmin, iinfo )
538  IF( iinfo.NE.0 ) THEN
539  info = 1
540  RETURN
541  END IF
542 *
543 * Choose Top-Down if D is (apparently) increasing,
544 * Bottom-Up if D is (apparently) decreasing.
545 *
546  IF( abs( d( 1 ) ).LE.abs( d( mnmin ) ) ) THEN
547  topdwn = .true.
548  ELSE
549  topdwn = .false.
550  END IF
551 *
552  IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
553 *
554 * Scale by DMAX
555 *
556  temp = abs( d( 1 ) )
557  DO 20 i = 2, mnmin
558  temp = max( temp, abs( d( i ) ) )
559  20 CONTINUE
560 *
561  IF( temp.GT.zero ) THEN
562  alpha = dmax / temp
563  ELSE
564  info = 2
565  RETURN
566  END IF
567 *
568  CALL sscal( mnmin, alpha, d, 1 )
569 *
570  END IF
571 *
572  CALL claset( 'Full', lda, n, czero, czero, a, lda )
573 *
574 * 3) Generate Banded Matrix using Givens rotations.
575 * Also the special case of UUB=LLB=0
576 *
577 * Compute Addressing constants to cover all
578 * storage formats. Whether GE, HE, SY, GB, HB, or SB,
579 * upper or lower triangle or both,
580 * the (i,j)-th element is in
581 * A( i - ISKEW*j + IOFFST, j )
582 *
583  IF( ipack.GT.4 ) THEN
584  ilda = lda - 1
585  iskew = 1
586  IF( ipack.GT.5 ) THEN
587  ioffst = uub + 1
588  ELSE
589  ioffst = 1
590  END IF
591  ELSE
592  ilda = lda
593  iskew = 0
594  ioffst = 0
595  END IF
596 *
597 * IPACKG is the format that the matrix is generated in. If this is
598 * different from IPACK, then the matrix must be repacked at the
599 * end. It also signals how to compute the norm, for scaling.
600 *
601  ipackg = 0
602 *
603 * Diagonal Matrix -- We are done, unless it
604 * is to be stored HP/SP/PP/TP (PACK='R' or 'C')
605 *
606  IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
607  DO 30 j = 1, mnmin
608  a( ( 1-iskew )*j+ioffst, j ) = cmplx( d( j ) )
609  30 CONTINUE
610 *
611  IF( ipack.LE.2 .OR. ipack.GE.5 )
612  $ ipackg = ipack
613 *
614  ELSE IF( givens ) THEN
615 *
616 * Check whether to use Givens rotations,
617 * Householder transformations, or nothing.
618 *
619  IF( isym.EQ.1 ) THEN
620 *
621 * Non-symmetric -- A = U D V
622 *
623  IF( ipack.GT.4 ) THEN
624  ipackg = ipack
625  ELSE
626  ipackg = 0
627  END IF
628 *
629  DO 40 j = 1, mnmin
630  a( ( 1-iskew )*j+ioffst, j ) = cmplx( d( j ) )
631  40 CONTINUE
632 *
633  IF( topdwn ) THEN
634  jkl = 0
635  DO 70 jku = 1, uub
636 *
637 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
638 *
639 * Last row actually rotated is M
640 * Last column actually rotated is MIN( M+JKU, N )
641 *
642  DO 60 jr = 1, min( m+jku, n ) + jkl - 1
643  extra = czero
644  angle = twopi*slarnd( 1, iseed )
645  c = cos( angle )*clarnd( 5, iseed )
646  s = sin( angle )*clarnd( 5, iseed )
647  icol = max( 1, jr-jkl )
648  IF( jr.LT.m ) THEN
649  il = min( n, jr+jku ) + 1 - icol
650  CALL clarot( .true., jr.GT.jkl, .false., il, c,
651  $ s, a( jr-iskew*icol+ioffst, icol ),
652  $ ilda, extra, dummy )
653  END IF
654 *
655 * Chase "EXTRA" back up
656 *
657  ir = jr
658  ic = icol
659  DO 50 jch = jr - jkl, 1, -jkl - jku
660  IF( ir.LT.m ) THEN
661  CALL clartg( a( ir+1-iskew*( ic+1 )+ioffst,
662  $ ic+1 ), extra, realc, s, dummy )
663  dummy = clarnd( 5, iseed )
664  c = conjg( realc*dummy )
665  s = conjg( -s*dummy )
666  END IF
667  irow = max( 1, jch-jku )
668  il = ir + 2 - irow
669  ctemp = czero
670  iltemp = jch.GT.jku
671  CALL clarot( .false., iltemp, .true., il, c, s,
672  $ a( irow-iskew*ic+ioffst, ic ),
673  $ ilda, ctemp, extra )
674  IF( iltemp ) THEN
675  CALL clartg( a( irow+1-iskew*( ic+1 )+ioffst,
676  $ ic+1 ), ctemp, realc, s, dummy )
677  dummy = clarnd( 5, iseed )
678  c = conjg( realc*dummy )
679  s = conjg( -s*dummy )
680 *
681  icol = max( 1, jch-jku-jkl )
682  il = ic + 2 - icol
683  extra = czero
684  CALL clarot( .true., jch.GT.jku+jkl, .true.,
685  $ il, c, s, a( irow-iskew*icol+
686  $ ioffst, icol ), ilda, extra,
687  $ ctemp )
688  ic = icol
689  ir = irow
690  END IF
691  50 CONTINUE
692  60 CONTINUE
693  70 CONTINUE
694 *
695  jku = uub
696  DO 100 jkl = 1, llb
697 *
698 * Transform from bandwidth JKL-1, JKU to JKL, JKU
699 *
700  DO 90 jc = 1, min( n+jkl, m ) + jku - 1
701  extra = czero
702  angle = twopi*slarnd( 1, iseed )
703  c = cos( angle )*clarnd( 5, iseed )
704  s = sin( angle )*clarnd( 5, iseed )
705  irow = max( 1, jc-jku )
706  IF( jc.LT.n ) THEN
707  il = min( m, jc+jkl ) + 1 - irow
708  CALL clarot( .false., jc.GT.jku, .false., il, c,
709  $ s, a( irow-iskew*jc+ioffst, jc ),
710  $ ilda, extra, dummy )
711  END IF
712 *
713 * Chase "EXTRA" back up
714 *
715  ic = jc
716  ir = irow
717  DO 80 jch = jc - jku, 1, -jkl - jku
718  IF( ic.LT.n ) THEN
719  CALL clartg( a( ir+1-iskew*( ic+1 )+ioffst,
720  $ ic+1 ), extra, realc, s, dummy )
721  dummy = clarnd( 5, iseed )
722  c = conjg( realc*dummy )
723  s = conjg( -s*dummy )
724  END IF
725  icol = max( 1, jch-jkl )
726  il = ic + 2 - icol
727  ctemp = czero
728  iltemp = jch.GT.jkl
729  CALL clarot( .true., iltemp, .true., il, c, s,
730  $ a( ir-iskew*icol+ioffst, icol ),
731  $ ilda, ctemp, extra )
732  IF( iltemp ) THEN
733  CALL clartg( a( ir+1-iskew*( icol+1 )+ioffst,
734  $ icol+1 ), ctemp, realc, s,
735  $ dummy )
736  dummy = clarnd( 5, iseed )
737  c = conjg( realc*dummy )
738  s = conjg( -s*dummy )
739  irow = max( 1, jch-jkl-jku )
740  il = ir + 2 - irow
741  extra = czero
742  CALL clarot( .false., jch.GT.jkl+jku, .true.,
743  $ il, c, s, a( irow-iskew*icol+
744  $ ioffst, icol ), ilda, extra,
745  $ ctemp )
746  ic = icol
747  ir = irow
748  END IF
749  80 CONTINUE
750  90 CONTINUE
751  100 CONTINUE
752 *
753  ELSE
754 *
755 * Bottom-Up -- Start at the bottom right.
756 *
757  jkl = 0
758  DO 130 jku = 1, uub
759 *
760 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
761 *
762 * First row actually rotated is M
763 * First column actually rotated is MIN( M+JKU, N )
764 *
765  iendch = min( m, n+jkl ) - 1
766  DO 120 jc = min( m+jku, n ) - 1, 1 - jkl, -1
767  extra = czero
768  angle = twopi*slarnd( 1, iseed )
769  c = cos( angle )*clarnd( 5, iseed )
770  s = sin( angle )*clarnd( 5, iseed )
771  irow = max( 1, jc-jku+1 )
772  IF( jc.GT.0 ) THEN
773  il = min( m, jc+jkl+1 ) + 1 - irow
774  CALL clarot( .false., .false., jc+jkl.LT.m, il,
775  $ c, s, a( irow-iskew*jc+ioffst,
776  $ jc ), ilda, dummy, extra )
777  END IF
778 *
779 * Chase "EXTRA" back down
780 *
781  ic = jc
782  DO 110 jch = jc + jkl, iendch, jkl + jku
783  ilextr = ic.GT.0
784  IF( ilextr ) THEN
785  CALL clartg( a( jch-iskew*ic+ioffst, ic ),
786  $ extra, realc, s, dummy )
787  dummy = clarnd( 5, iseed )
788  c = realc*dummy
789  s = s*dummy
790  END IF
791  ic = max( 1, ic )
792  icol = min( n-1, jch+jku )
793  iltemp = jch + jku.LT.n
794  ctemp = czero
795  CALL clarot( .true., ilextr, iltemp, icol+2-ic,
796  $ c, s, a( jch-iskew*ic+ioffst, ic ),
797  $ ilda, extra, ctemp )
798  IF( iltemp ) THEN
799  CALL clartg( a( jch-iskew*icol+ioffst,
800  $ icol ), ctemp, realc, s, dummy )
801  dummy = clarnd( 5, iseed )
802  c = realc*dummy
803  s = s*dummy
804  il = min( iendch, jch+jkl+jku ) + 2 - jch
805  extra = czero
806  CALL clarot( .false., .true.,
807  $ jch+jkl+jku.LE.iendch, il, c, s,
808  $ a( jch-iskew*icol+ioffst,
809  $ icol ), ilda, ctemp, extra )
810  ic = icol
811  END IF
812  110 CONTINUE
813  120 CONTINUE
814  130 CONTINUE
815 *
816  jku = uub
817  DO 160 jkl = 1, llb
818 *
819 * Transform from bandwidth JKL-1, JKU to JKL, JKU
820 *
821 * First row actually rotated is MIN( N+JKL, M )
822 * First column actually rotated is N
823 *
824  iendch = min( n, m+jku ) - 1
825  DO 150 jr = min( n+jkl, m ) - 1, 1 - jku, -1
826  extra = czero
827  angle = twopi*slarnd( 1, iseed )
828  c = cos( angle )*clarnd( 5, iseed )
829  s = sin( angle )*clarnd( 5, iseed )
830  icol = max( 1, jr-jkl+1 )
831  IF( jr.GT.0 ) THEN
832  il = min( n, jr+jku+1 ) + 1 - icol
833  CALL clarot( .true., .false., jr+jku.LT.n, il,
834  $ c, s, a( jr-iskew*icol+ioffst,
835  $ icol ), ilda, dummy, extra )
836  END IF
837 *
838 * Chase "EXTRA" back down
839 *
840  ir = jr
841  DO 140 jch = jr + jku, iendch, jkl + jku
842  ilextr = ir.GT.0
843  IF( ilextr ) THEN
844  CALL clartg( a( ir-iskew*jch+ioffst, jch ),
845  $ extra, realc, s, dummy )
846  dummy = clarnd( 5, iseed )
847  c = realc*dummy
848  s = s*dummy
849  END IF
850  ir = max( 1, ir )
851  irow = min( m-1, jch+jkl )
852  iltemp = jch + jkl.LT.m
853  ctemp = czero
854  CALL clarot( .false., ilextr, iltemp, irow+2-ir,
855  $ c, s, a( ir-iskew*jch+ioffst,
856  $ jch ), ilda, extra, ctemp )
857  IF( iltemp ) THEN
858  CALL clartg( a( irow-iskew*jch+ioffst, jch ),
859  $ ctemp, realc, s, dummy )
860  dummy = clarnd( 5, iseed )
861  c = realc*dummy
862  s = s*dummy
863  il = min( iendch, jch+jkl+jku ) + 2 - jch
864  extra = czero
865  CALL clarot( .true., .true.,
866  $ jch+jkl+jku.LE.iendch, il, c, s,
867  $ a( irow-iskew*jch+ioffst, jch ),
868  $ ilda, ctemp, extra )
869  ir = irow
870  END IF
871  140 CONTINUE
872  150 CONTINUE
873  160 CONTINUE
874 *
875  END IF
876 *
877  ELSE
878 *
879 * Symmetric -- A = U D U'
880 * Hermitian -- A = U D U*
881 *
882  ipackg = ipack
883  ioffg = ioffst
884 *
885  IF( topdwn ) THEN
886 *
887 * Top-Down -- Generate Upper triangle only
888 *
889  IF( ipack.GE.5 ) THEN
890  ipackg = 6
891  ioffg = uub + 1
892  ELSE
893  ipackg = 1
894  END IF
895 *
896  DO 170 j = 1, mnmin
897  a( ( 1-iskew )*j+ioffg, j ) = cmplx( d( j ) )
898  170 CONTINUE
899 *
900  DO 200 k = 1, uub
901  DO 190 jc = 1, n - 1
902  irow = max( 1, jc-k )
903  il = min( jc+1, k+2 )
904  extra = czero
905  ctemp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
906  angle = twopi*slarnd( 1, iseed )
907  c = cos( angle )*clarnd( 5, iseed )
908  s = sin( angle )*clarnd( 5, iseed )
909  IF( csym ) THEN
910  ct = c
911  st = s
912  ELSE
913  ctemp = conjg( ctemp )
914  ct = conjg( c )
915  st = conjg( s )
916  END IF
917  CALL clarot( .false., jc.GT.k, .true., il, c, s,
918  $ a( irow-iskew*jc+ioffg, jc ), ilda,
919  $ extra, ctemp )
920  CALL clarot( .true., .true., .false.,
921  $ min( k, n-jc )+1, ct, st,
922  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
923  $ ctemp, dummy )
924 *
925 * Chase EXTRA back up the matrix
926 *
927  icol = jc
928  DO 180 jch = jc - k, 1, -k
929  CALL clartg( a( jch+1-iskew*( icol+1 )+ioffg,
930  $ icol+1 ), extra, realc, s, dummy )
931  dummy = clarnd( 5, iseed )
932  c = conjg( realc*dummy )
933  s = conjg( -s*dummy )
934  ctemp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
935  IF( csym ) THEN
936  ct = c
937  st = s
938  ELSE
939  ctemp = conjg( ctemp )
940  ct = conjg( c )
941  st = conjg( s )
942  END IF
943  CALL clarot( .true., .true., .true., k+2, c, s,
944  $ a( ( 1-iskew )*jch+ioffg, jch ),
945  $ ilda, ctemp, extra )
946  irow = max( 1, jch-k )
947  il = min( jch+1, k+2 )
948  extra = czero
949  CALL clarot( .false., jch.GT.k, .true., il, ct,
950  $ st, a( irow-iskew*jch+ioffg, jch ),
951  $ ilda, extra, ctemp )
952  icol = jch
953  180 CONTINUE
954  190 CONTINUE
955  200 CONTINUE
956 *
957 * If we need lower triangle, copy from upper. Note that
958 * the order of copying is chosen to work for 'q' -> 'b'
959 *
960  IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
961  DO 230 jc = 1, n
962  irow = ioffst - iskew*jc
963  IF( csym ) THEN
964  DO 210 jr = jc, min( n, jc+uub )
965  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
966  210 CONTINUE
967  ELSE
968  DO 220 jr = jc, min( n, jc+uub )
969  a( jr+irow, jc ) = conjg( a( jc-iskew*jr+
970  $ ioffg, jr ) )
971  220 CONTINUE
972  END IF
973  230 CONTINUE
974  IF( ipack.EQ.5 ) THEN
975  DO 250 jc = n - uub + 1, n
976  DO 240 jr = n + 2 - jc, uub + 1
977  a( jr, jc ) = czero
978  240 CONTINUE
979  250 CONTINUE
980  END IF
981  IF( ipackg.EQ.6 ) THEN
982  ipackg = ipack
983  ELSE
984  ipackg = 0
985  END IF
986  END IF
987  ELSE
988 *
989 * Bottom-Up -- Generate Lower triangle only
990 *
991  IF( ipack.GE.5 ) THEN
992  ipackg = 5
993  IF( ipack.EQ.6 )
994  $ ioffg = 1
995  ELSE
996  ipackg = 2
997  END IF
998 *
999  DO 260 j = 1, mnmin
1000  a( ( 1-iskew )*j+ioffg, j ) = cmplx( d( j ) )
1001  260 CONTINUE
1002 *
1003  DO 290 k = 1, uub
1004  DO 280 jc = n - 1, 1, -1
1005  il = min( n+1-jc, k+2 )
1006  extra = czero
1007  ctemp = a( 1+( 1-iskew )*jc+ioffg, jc )
1008  angle = twopi*slarnd( 1, iseed )
1009  c = cos( angle )*clarnd( 5, iseed )
1010  s = sin( angle )*clarnd( 5, iseed )
1011  IF( csym ) THEN
1012  ct = c
1013  st = s
1014  ELSE
1015  ctemp = conjg( ctemp )
1016  ct = conjg( c )
1017  st = conjg( s )
1018  END IF
1019  CALL clarot( .false., .true., n-jc.GT.k, il, c, s,
1020  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
1021  $ ctemp, extra )
1022  icol = max( 1, jc-k+1 )
1023  CALL clarot( .true., .false., .true., jc+2-icol,
1024  $ ct, st, a( jc-iskew*icol+ioffg,
1025  $ icol ), ilda, dummy, ctemp )
1026 *
1027 * Chase EXTRA back down the matrix
1028 *
1029  icol = jc
1030  DO 270 jch = jc + k, n - 1, k
1031  CALL clartg( a( jch-iskew*icol+ioffg, icol ),
1032  $ extra, realc, s, dummy )
1033  dummy = clarnd( 5, iseed )
1034  c = realc*dummy
1035  s = s*dummy
1036  ctemp = a( 1+( 1-iskew )*jch+ioffg, jch )
1037  IF( csym ) THEN
1038  ct = c
1039  st = s
1040  ELSE
1041  ctemp = conjg( ctemp )
1042  ct = conjg( c )
1043  st = conjg( s )
1044  END IF
1045  CALL clarot( .true., .true., .true., k+2, c, s,
1046  $ a( jch-iskew*icol+ioffg, icol ),
1047  $ ilda, extra, ctemp )
1048  il = min( n+1-jch, k+2 )
1049  extra = czero
1050  CALL clarot( .false., .true., n-jch.GT.k, il,
1051  $ ct, st, a( ( 1-iskew )*jch+ioffg,
1052  $ jch ), ilda, ctemp, extra )
1053  icol = jch
1054  270 CONTINUE
1055  280 CONTINUE
1056  290 CONTINUE
1057 *
1058 * If we need upper triangle, copy from lower. Note that
1059 * the order of copying is chosen to work for 'b' -> 'q'
1060 *
1061  IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
1062  DO 320 jc = n, 1, -1
1063  irow = ioffst - iskew*jc
1064  IF( csym ) THEN
1065  DO 300 jr = jc, max( 1, jc-uub ), -1
1066  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
1067  300 CONTINUE
1068  ELSE
1069  DO 310 jr = jc, max( 1, jc-uub ), -1
1070  a( jr+irow, jc ) = conjg( a( jc-iskew*jr+
1071  $ ioffg, jr ) )
1072  310 CONTINUE
1073  END IF
1074  320 CONTINUE
1075  IF( ipack.EQ.6 ) THEN
1076  DO 340 jc = 1, uub
1077  DO 330 jr = 1, uub + 1 - jc
1078  a( jr, jc ) = czero
1079  330 CONTINUE
1080  340 CONTINUE
1081  END IF
1082  IF( ipackg.EQ.5 ) THEN
1083  ipackg = ipack
1084  ELSE
1085  ipackg = 0
1086  END IF
1087  END IF
1088  END IF
1089 *
1090 * Ensure that the diagonal is real if Hermitian
1091 *
1092  IF( .NOT.csym ) THEN
1093  DO 350 jc = 1, n
1094  irow = ioffst + ( 1-iskew )*jc
1095  a( irow, jc ) = cmplx( real( a( irow, jc ) ) )
1096  350 CONTINUE
1097  END IF
1098 *
1099  END IF
1100 *
1101  ELSE
1102 *
1103 * 4) Generate Banded Matrix by first
1104 * Rotating by random Unitary matrices,
1105 * then reducing the bandwidth using Householder
1106 * transformations.
1107 *
1108 * Note: we should get here only if LDA .ge. N
1109 *
1110  IF( isym.EQ.1 ) THEN
1111 *
1112 * Non-symmetric -- A = U D V
1113 *
1114  CALL clagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1115  $ iinfo )
1116  ELSE
1117 *
1118 * Symmetric -- A = U D U' or
1119 * Hermitian -- A = U D U*
1120 *
1121  IF( csym ) THEN
1122  CALL clagsy( m, llb, d, a, lda, iseed, work, iinfo )
1123  ELSE
1124  CALL claghe( m, llb, d, a, lda, iseed, work, iinfo )
1125  END IF
1126  END IF
1127 *
1128  IF( iinfo.NE.0 ) THEN
1129  info = 3
1130  RETURN
1131  END IF
1132  END IF
1133 *
1134 * 5) Pack the matrix
1135 *
1136  IF( ipack.NE.ipackg ) THEN
1137  IF( ipack.EQ.1 ) THEN
1138 *
1139 * 'U' -- Upper triangular, not packed
1140 *
1141  DO 370 j = 1, m
1142  DO 360 i = j + 1, m
1143  a( i, j ) = czero
1144  360 CONTINUE
1145  370 CONTINUE
1146 *
1147  ELSE IF( ipack.EQ.2 ) THEN
1148 *
1149 * 'L' -- Lower triangular, not packed
1150 *
1151  DO 390 j = 2, m
1152  DO 380 i = 1, j - 1
1153  a( i, j ) = czero
1154  380 CONTINUE
1155  390 CONTINUE
1156 *
1157  ELSE IF( ipack.EQ.3 ) THEN
1158 *
1159 * 'C' -- Upper triangle packed Columnwise.
1160 *
1161  icol = 1
1162  irow = 0
1163  DO 410 j = 1, m
1164  DO 400 i = 1, j
1165  irow = irow + 1
1166  IF( irow.GT.lda ) THEN
1167  irow = 1
1168  icol = icol + 1
1169  END IF
1170  a( irow, icol ) = a( i, j )
1171  400 CONTINUE
1172  410 CONTINUE
1173 *
1174  ELSE IF( ipack.EQ.4 ) THEN
1175 *
1176 * 'R' -- Lower triangle packed Columnwise.
1177 *
1178  icol = 1
1179  irow = 0
1180  DO 430 j = 1, m
1181  DO 420 i = j, m
1182  irow = irow + 1
1183  IF( irow.GT.lda ) THEN
1184  irow = 1
1185  icol = icol + 1
1186  END IF
1187  a( irow, icol ) = a( i, j )
1188  420 CONTINUE
1189  430 CONTINUE
1190 *
1191  ELSE IF( ipack.GE.5 ) THEN
1192 *
1193 * 'B' -- The lower triangle is packed as a band matrix.
1194 * 'Q' -- The upper triangle is packed as a band matrix.
1195 * 'Z' -- The whole matrix is packed as a band matrix.
1196 *
1197  IF( ipack.EQ.5 )
1198  $ uub = 0
1199  IF( ipack.EQ.6 )
1200  $ llb = 0
1201 *
1202  DO 450 j = 1, uub
1203  DO 440 i = min( j+llb, m ), 1, -1
1204  a( i-j+uub+1, j ) = a( i, j )
1205  440 CONTINUE
1206  450 CONTINUE
1207 *
1208  DO 470 j = uub + 2, n
1209  DO 460 i = j - uub, min( j+llb, m )
1210  a( i-j+uub+1, j ) = a( i, j )
1211  460 CONTINUE
1212  470 CONTINUE
1213  END IF
1214 *
1215 * If packed, zero out extraneous elements.
1216 *
1217 * Symmetric/Triangular Packed --
1218 * zero out everything after A(IROW,ICOL)
1219 *
1220  IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1221  DO 490 jc = icol, m
1222  DO 480 jr = irow + 1, lda
1223  a( jr, jc ) = czero
1224  480 CONTINUE
1225  irow = 0
1226  490 CONTINUE
1227 *
1228  ELSE IF( ipack.GE.5 ) THEN
1229 *
1230 * Packed Band --
1231 * 1st row is now in A( UUB+2-j, j), zero above it
1232 * m-th row is now in A( M+UUB-j,j), zero below it
1233 * last non-zero diagonal is now in A( UUB+LLB+1,j ),
1234 * zero below it, too.
1235 *
1236  ir1 = uub + llb + 2
1237  ir2 = uub + m + 2
1238  DO 520 jc = 1, n
1239  DO 500 jr = 1, uub + 1 - jc
1240  a( jr, jc ) = czero
1241  500 CONTINUE
1242  DO 510 jr = max( 1, min( ir1, ir2-jc ) ), lda
1243  a( jr, jc ) = czero
1244  510 CONTINUE
1245  520 CONTINUE
1246  END IF
1247  END IF
1248 *
1249  RETURN
1250 *
1251 * End of CLATMS
1252 *
1253  END
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine clagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
CLAGGE
Definition: clagge.f:114
subroutine clatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
CLATMS
Definition: clatms.f:332
subroutine claghe(N, K, D, A, LDA, ISEED, WORK, INFO)
CLAGHE
Definition: claghe.f:102
subroutine clarot(LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
CLAROT
Definition: clarot.f:229
subroutine clagsy(N, K, D, A, LDA, ISEED, WORK, INFO)
CLAGSY
Definition: clagsy.f:102
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine slatm1(MODE, COND, IRSIGN, IDIST, ISEED, D, N, INFO)
SLATM1
Definition: slatm1.f:135
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79