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