SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slatms()

subroutine slatms ( integer  m,
integer  n,
character  dist,
integer, dimension( 4 )  iseed,
character  sym,
real, dimension( * )  d,
integer  mode,
real  cond,
real  dmax,
integer  kl,
integer  ku,
character  pack,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( * )  work,
integer  info 
)

Definition at line 1 of file slatms.f.

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