LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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,
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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine clagge(m, n, kl, ku, d, a, lda, iseed, work, info)
CLAGGE
Definition clagge.f:114
subroutine claghe(n, k, d, a, lda, iseed, work, info)
CLAGHE
Definition claghe.f:102
subroutine clagsy(n, k, d, a, lda, iseed, work, info)
CLAGSY
Definition clagsy.f:102
subroutine clarot(lrows, lleft, lright, nl, c, s, a, lda, xleft, xright)
CLAROT
Definition clarot.f:229
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
Definition clatms.f:332
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:116
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 sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine slatm1(mode, cond, irsign, idist, iseed, d, n, info)
SLATM1
Definition slatm1.f:135