LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slatmt()

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

SLATMT

Purpose:
!>
!>    SLATMT generates random matrices with specified singular values
!>    (or symmetric/hermitian with specified eigenvalues)
!>    for testing LAPACK programs.
!>
!>    SLATMT operates by applying the following sequence of
!>    operations:
!>
!>      Set the diagonal to D, where D may be input or
!>         computed according to MODE, COND, DMAX, and SYM
!>         as described below.
!>
!>      Generate a matrix with the appropriate band structure, by one
!>         of two methods:
!>
!>      Method A:
!>          Generate a dense M x N matrix by multiplying D on the left
!>              and the right by random unitary matrices, then:
!>
!>          Reduce the bandwidth according to KL and KU, using
!>          Householder transformations.
!>
!>      Method B:
!>          Convert the bandwidth-0 (i.e., diagonal) matrix to a
!>              bandwidth-1 matrix using Givens rotations, 
!>              out-of-band elements back, much as in QR; then
!>              convert the bandwidth-1 to a bandwidth-2 matrix, etc.
!>              Note that for reasonably small bandwidths (relative to
!>              M and N) this requires less storage, as a dense matrix
!>              is not generated.  Also, for symmetric matrices, only
!>              one triangle is generated.
!>
!>      Method A is chosen if the bandwidth is a large fraction of the
!>          order of the matrix, and LDA is at least M (so a dense
!>          matrix can be stored.)  Method B is chosen if the bandwidth
!>          is small (< 1/2 N for symmetric, < .3 N+M for
!>          non-symmetric), or LDA is less than M and not less than the
!>          bandwidth.
!>
!>      Pack the matrix if desired. Options specified by PACK are:
!>         no packing
!>         zero out upper half (if symmetric)
!>         zero out lower half (if symmetric)
!>         store the upper half columnwise (if symmetric or upper
!>               triangular)
!>         store the lower half columnwise (if symmetric or lower
!>               triangular)
!>         store the lower triangle in banded format (if symmetric
!>               or lower triangular)
!>         store the upper triangle in banded format (if symmetric
!>               or upper triangular)
!>         store the entire matrix in banded format
!>      If Method B is chosen, and band format is specified, then the
!>         matrix will be generated in the band format, so no repacking
!>         will be necessary.
!> 
Parameters
[in]M
!>          M is INTEGER
!>           The number of rows of A. Not modified.
!> 
[in]N
!>          N is INTEGER
!>           The number of columns of A. Not modified.
!> 
[in]DIST
!>          DIST is CHARACTER*1
!>           On entry, DIST specifies the type of distribution to be used
!>           to generate the random eigen-/singular values.
!>           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
!>           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
!>           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
!>           Not modified.
!> 
[in,out]ISEED
!>          ISEED is INTEGER array, dimension ( 4 )
!>           On entry ISEED specifies the seed of the random number
!>           generator. They should lie between 0 and 4095 inclusive,
!>           and ISEED(4) should be odd. The random number generator
!>           uses a linear congruential sequence limited to small
!>           integers, and so should produce machine independent
!>           random numbers. The values of ISEED are changed on
!>           exit, and can be used in the next call to SLATMT
!>           to continue the same random number sequence.
!>           Changed on exit.
!> 
[in]SYM
!>          SYM is CHARACTER*1
!>           If SYM='S' or 'H', the generated matrix is symmetric, with
!>             eigenvalues specified by D, COND, MODE, and DMAX; they
!>             may be positive, negative, or zero.
!>           If SYM='P', the generated matrix is symmetric, with
!>             eigenvalues (= singular values) specified by D, COND,
!>             MODE, and DMAX; they will not be negative.
!>           If SYM='N', the generated matrix is nonsymmetric, with
!>             singular values specified by D, COND, MODE, and DMAX;
!>             they will not be negative.
!>           Not modified.
!> 
[in,out]D
!>          D is REAL array, dimension ( MIN( M , N ) )
!>           This array is used to specify the singular values or
!>           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
!>           assumed to contain the singular/eigenvalues, otherwise
!>           they will be computed according to MODE, COND, and DMAX,
!>           and placed in D.
!>           Modified if MODE is nonzero.
!> 
[in]MODE
!>          MODE is INTEGER
!>           On entry this describes how the singular/eigenvalues are to
!>           be specified:
!>           MODE = 0 means use D as input
!>
!>           MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
!>           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
!>           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
!>
!>           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
!>           MODE = 5 sets D to random numbers in the range
!>                    ( 1/COND , 1 ) such that their logarithms
!>                    are uniformly distributed.
!>           MODE = 6 set D to random numbers from same distribution
!>                    as the rest of the matrix.
!>           MODE < 0 has the same meaning as ABS(MODE), except that
!>              the order of the elements of D is reversed.
!>           Thus if MODE is positive, D has entries ranging from
!>              1 to 1/COND, if negative, from 1/COND to 1,
!>           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
!>              the elements of D will also be multiplied by a random
!>              sign (i.e., +1 or -1.)
!>           Not modified.
!> 
[in]COND
!>          COND is REAL
!>           On entry, this is used as described under MODE above.
!>           If used, it must be >= 1. Not modified.
!> 
[in]DMAX
!>          DMAX is REAL
!>           If MODE is neither -6, 0 nor 6, the contents of D, as
!>           computed according to MODE and COND, will be scaled by
!>           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
!>           singular value (which is to say the norm) will be abs(DMAX).
!>           Note that DMAX need not be positive: if DMAX is negative
!>           (or zero), D will be scaled by a negative number (or zero).
!>           Not modified.
!> 
[in]RANK
!>          RANK is INTEGER
!>           The rank of matrix to be generated for modes 1,2,3 only.
!>           D( RANK+1:N ) = 0.
!>           Not modified.
!> 
[in]KL
!>          KL is INTEGER
!>           This specifies the lower bandwidth of the  matrix. For
!>           example, KL=0 implies upper triangular, KL=1 implies upper
!>           Hessenberg, and KL being at least M-1 means that the matrix
!>           has full lower bandwidth.  KL must equal KU if the matrix
!>           is symmetric.
!>           Not modified.
!> 
[in]KU
!>          KU is INTEGER
!>           This specifies the upper bandwidth of the  matrix. For
!>           example, KU=0 implies lower triangular, KU=1 implies lower
!>           Hessenberg, and KU being at least N-1 means that the matrix
!>           has full upper bandwidth.  KL must equal KU if the matrix
!>           is symmetric.
!>           Not modified.
!> 
[in]PACK
!>          PACK is CHARACTER*1
!>           This specifies packing of matrix as follows:
!>           'N' => no packing
!>           'U' => zero out all subdiagonal entries (if symmetric)
!>           'L' => zero out all superdiagonal entries (if symmetric)
!>           'C' => store the upper triangle columnwise
!>                  (only if the matrix is symmetric or upper triangular)
!>           'R' => store the lower triangle columnwise
!>                  (only if the matrix is symmetric or lower triangular)
!>           'B' => store the lower triangle in band storage scheme
!>                  (only if matrix symmetric or lower triangular)
!>           'Q' => store the upper triangle in band storage scheme
!>                  (only if matrix symmetric or upper triangular)
!>           'Z' => store the entire matrix in band storage scheme
!>                      (pivoting can be provided for by using this
!>                      option to store A in the trailing rows of
!>                      the allocated storage)
!>
!>           Using these options, the various LAPACK packed and banded
!>           storage schemes can be obtained:
!>           GB               - use 'Z'
!>           PB, SB or TB     - use 'B' or 'Q'
!>           PP, SP or TP     - use 'C' or 'R'
!>
!>           If two calls to SLATMT differ only in the PACK parameter,
!>           they will generate mathematically equivalent matrices.
!>           Not modified.
!> 
[in,out]A
!>          A is REAL array, dimension ( LDA, N )
!>           On exit A is the desired test matrix.  A is first generated
!>           in full (unpacked) form, and then packed, if so specified
!>           by PACK.  Thus, the first M elements of the first N
!>           columns will always be modified.  If PACK specifies a
!>           packed or banded storage scheme, all LDA elements of the
!>           first N columns will be modified; the elements of the
!>           array which do not correspond to elements of the generated
!>           matrix are set to zero.
!>           Modified.
!> 
[in]LDA
!>          LDA is INTEGER
!>           LDA specifies the first dimension of A as declared in the
!>           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
!>           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
!>           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
!>           If PACK='Z', LDA must be large enough to hold the packed
!>           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
!>           Not modified.
!> 
[out]WORK
!>          WORK is REAL array, dimension ( 3*MAX( N , M ) )
!>           Workspace.
!>           Modified.
!> 
[out]INFO
!>          INFO is INTEGER
!>           Error code.  On exit, INFO will be set to one of the
!>           following values:
!>             0 => normal return
!>            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
!>            -2 => N negative
!>            -3 => DIST illegal string
!>            -5 => SYM illegal string
!>            -7 => MODE not in range -6 to 6
!>            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
!>           -10 => KL negative
!>           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
!>           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
!>                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
!>                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
!>                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
!>                  N.
!>           -14 => LDA is less than M, or PACK='Z' and LDA is less than
!>                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
!>            1  => Error return from SLATM7
!>            2  => Cannot scale to DMAX (max. sing. value is 0)
!>            3  => Error return from SLAGGE or SLAGSY
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 329 of file slatmt.f.

331*
332* -- LAPACK computational routine --
333* -- LAPACK is a software package provided by Univ. of Tennessee, --
334* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335*
336* .. Scalar Arguments ..
337 REAL COND, DMAX
338 INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK
339 CHARACTER DIST, PACK, SYM
340* ..
341* .. Array Arguments ..
342 REAL A( LDA, * ), D( * ), WORK( * )
343 INTEGER ISEED( 4 )
344* ..
345*
346* =====================================================================
347*
348* .. Parameters ..
349 REAL ZERO
350 parameter( zero = 0.0e0 )
351 REAL ONE
352 parameter( one = 1.0e0 )
353 REAL TWOPI
354 parameter( twopi = 6.28318530717958647692528676655900576839e+0 )
355* ..
356* .. Local Scalars ..
357 REAL ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
358 INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
359 $ IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2,
360 $ IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH,
361 $ JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC,
362 $ UUB
363 LOGICAL GIVENS, ILEXTR, ILTEMP, TOPDWN
364* ..
365* .. External Functions ..
366 REAL SLARND
367 LOGICAL LSAME
368 EXTERNAL slarnd, lsame
369* ..
370* .. External Subroutines ..
371 EXTERNAL slatm7, scopy, slagge, slagsy, slarot,
373* ..
374* .. Intrinsic Functions ..
375 INTRINSIC abs, cos, max, min, mod, real, sin
376* ..
377* .. Executable Statements ..
378*
379* 1) Decode and Test the input parameters.
380* Initialize flags & seed.
381*
382 info = 0
383*
384* Quick return if possible
385*
386 IF( m.EQ.0 .OR. n.EQ.0 )
387 $ RETURN
388*
389* Decode DIST
390*
391 IF( lsame( dist, 'U' ) ) THEN
392 idist = 1
393 ELSE IF( lsame( dist, 'S' ) ) THEN
394 idist = 2
395 ELSE IF( lsame( dist, 'N' ) ) THEN
396 idist = 3
397 ELSE
398 idist = -1
399 END IF
400*
401* Decode SYM
402*
403 IF( lsame( sym, 'N' ) ) THEN
404 isym = 1
405 irsign = 0
406 ELSE IF( lsame( sym, 'P' ) ) THEN
407 isym = 2
408 irsign = 0
409 ELSE IF( lsame( sym, 'S' ) ) THEN
410 isym = 2
411 irsign = 1
412 ELSE IF( lsame( sym, 'H' ) ) THEN
413 isym = 2
414 irsign = 1
415 ELSE
416 isym = -1
417 END IF
418*
419* Decode PACK
420*
421 isympk = 0
422 IF( lsame( pack, 'N' ) ) THEN
423 ipack = 0
424 ELSE IF( lsame( pack, 'U' ) ) THEN
425 ipack = 1
426 isympk = 1
427 ELSE IF( lsame( pack, 'L' ) ) THEN
428 ipack = 2
429 isympk = 1
430 ELSE IF( lsame( pack, 'C' ) ) THEN
431 ipack = 3
432 isympk = 2
433 ELSE IF( lsame( pack, 'R' ) ) THEN
434 ipack = 4
435 isympk = 3
436 ELSE IF( lsame( pack, 'B' ) ) THEN
437 ipack = 5
438 isympk = 3
439 ELSE IF( lsame( pack, 'Q' ) ) THEN
440 ipack = 6
441 isympk = 2
442 ELSE IF( lsame( pack, 'Z' ) ) THEN
443 ipack = 7
444 ELSE
445 ipack = -1
446 END IF
447*
448* Set certain internal parameters
449*
450 mnmin = min( m, n )
451 llb = min( kl, m-1 )
452 uub = min( ku, n-1 )
453 mr = min( m, n+llb )
454 nc = min( n, m+uub )
455*
456 IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
457 minlda = uub + 1
458 ELSE IF( ipack.EQ.7 ) THEN
459 minlda = llb + uub + 1
460 ELSE
461 minlda = m
462 END IF
463*
464* Use Givens rotation method if bandwidth small enough,
465* or if LDA is too small to store the matrix unpacked.
466*
467 givens = .false.
468 IF( isym.EQ.1 ) THEN
469 IF( real( llb+uub ).LT.0.3*real( max( 1, mr+nc ) ) )
470 $ givens = .true.
471 ELSE
472 IF( 2*llb.LT.m )
473 $ givens = .true.
474 END IF
475 IF( lda.LT.m .AND. lda.GE.minlda )
476 $ givens = .true.
477*
478* Set INFO if an error
479*
480 IF( m.LT.0 ) THEN
481 info = -1
482 ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
483 info = -1
484 ELSE IF( n.LT.0 ) THEN
485 info = -2
486 ELSE IF( idist.EQ.-1 ) THEN
487 info = -3
488 ELSE IF( isym.EQ.-1 ) THEN
489 info = -5
490 ELSE IF( abs( mode ).GT.6 ) THEN
491 info = -7
492 ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
493 $ THEN
494 info = -8
495 ELSE IF( kl.LT.0 ) THEN
496 info = -10
497 ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
498 info = -11
499 ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
500 $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
501 $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
502 $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
503 info = -12
504 ELSE IF( lda.LT.max( 1, minlda ) ) THEN
505 info = -14
506 END IF
507*
508 IF( info.NE.0 ) THEN
509 CALL xerbla( 'SLATMT', -info )
510 RETURN
511 END IF
512*
513* Initialize random number generator
514*
515 DO 100 i = 1, 4
516 iseed( i ) = mod( abs( iseed( i ) ), 4096 )
517 100 CONTINUE
518*
519 IF( mod( iseed( 4 ), 2 ).NE.1 )
520 $ iseed( 4 ) = iseed( 4 ) + 1
521*
522* 2) Set up D if indicated.
523*
524* Compute D according to COND and MODE
525*
526 CALL slatm7( mode, cond, irsign, idist, iseed, d, mnmin, rank,
527 $ iinfo )
528 IF( iinfo.NE.0 ) THEN
529 info = 1
530 RETURN
531 END IF
532*
533* Choose Top-Down if D is (apparently) increasing,
534* Bottom-Up if D is (apparently) decreasing.
535*
536 IF( abs( d( 1 ) ).LE.abs( d( rank ) ) ) THEN
537 topdwn = .true.
538 ELSE
539 topdwn = .false.
540 END IF
541*
542 IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
543*
544* Scale by DMAX
545*
546 temp = abs( d( 1 ) )
547 DO 110 i = 2, rank
548 temp = max( temp, abs( d( i ) ) )
549 110 CONTINUE
550*
551 IF( temp.GT.zero ) THEN
552 alpha = dmax / temp
553 ELSE
554 info = 2
555 RETURN
556 END IF
557*
558 CALL sscal( rank, alpha, d, 1 )
559*
560 END IF
561*
562* 3) Generate Banded Matrix using Givens rotations.
563* Also the special case of UUB=LLB=0
564*
565* Compute Addressing constants to cover all
566* storage formats. Whether GE, SY, GB, or SB,
567* upper or lower triangle or both,
568* the (i,j)-th element is in
569* A( i - ISKEW*j + IOFFST, j )
570*
571 IF( ipack.GT.4 ) THEN
572 ilda = lda - 1
573 iskew = 1
574 IF( ipack.GT.5 ) THEN
575 ioffst = uub + 1
576 ELSE
577 ioffst = 1
578 END IF
579 ELSE
580 ilda = lda
581 iskew = 0
582 ioffst = 0
583 END IF
584*
585* IPACKG is the format that the matrix is generated in. If this is
586* different from IPACK, then the matrix must be repacked at the
587* end. It also signals how to compute the norm, for scaling.
588*
589 ipackg = 0
590 CALL slaset( 'Full', lda, n, zero, zero, a, lda )
591*
592* Diagonal Matrix -- We are done, unless it
593* is to be stored SP/PP/TP (PACK='R' or 'C')
594*
595 IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
596 CALL scopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
597 IF( ipack.LE.2 .OR. ipack.GE.5 )
598 $ ipackg = ipack
599*
600 ELSE IF( givens ) THEN
601*
602* Check whether to use Givens rotations,
603* Householder transformations, or nothing.
604*
605 IF( isym.EQ.1 ) THEN
606*
607* Non-symmetric -- A = U D V
608*
609 IF( ipack.GT.4 ) THEN
610 ipackg = ipack
611 ELSE
612 ipackg = 0
613 END IF
614*
615 CALL scopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
616*
617 IF( topdwn ) THEN
618 jkl = 0
619 DO 140 jku = 1, uub
620*
621* Transform from bandwidth JKL, JKU-1 to JKL, JKU
622*
623* Last row actually rotated is M
624* Last column actually rotated is MIN( M+JKU, N )
625*
626 DO 130 jr = 1, min( m+jku, n ) + jkl - 1
627 extra = zero
628 angle = twopi*slarnd( 1, iseed )
629 c = cos( angle )
630 s = sin( angle )
631 icol = max( 1, jr-jkl )
632 IF( jr.LT.m ) THEN
633 il = min( n, jr+jku ) + 1 - icol
634 CALL slarot( .true., jr.GT.jkl, .false., il,
635 $ c,
636 $ s, a( jr-iskew*icol+ioffst, icol ),
637 $ ilda, extra, dummy )
638 END IF
639*
640* Chase "EXTRA" back up
641*
642 ir = jr
643 ic = icol
644 DO 120 jch = jr - jkl, 1, -jkl - jku
645 IF( ir.LT.m ) THEN
646 CALL slartg( a( ir+1-iskew*( ic+1 )+ioffst,
647 $ ic+1 ), extra, c, s, dummy )
648 END IF
649 irow = max( 1, jch-jku )
650 il = ir + 2 - irow
651 temp = zero
652 iltemp = jch.GT.jku
653 CALL slarot( .false., iltemp, .true., il, c,
654 $ -s,
655 $ a( irow-iskew*ic+ioffst, ic ),
656 $ ilda, temp, extra )
657 IF( iltemp ) THEN
658 CALL slartg( a( irow+1-iskew*( ic+1 )+ioffst,
659 $ ic+1 ), temp, c, s, dummy )
660 icol = max( 1, jch-jku-jkl )
661 il = ic + 2 - icol
662 extra = zero
663 CALL slarot( .true., jch.GT.jku+jkl,
664 $ .true.,
665 $ il, c, -s, a( irow-iskew*icol+
666 $ ioffst, icol ), ilda, extra,
667 $ temp )
668 ic = icol
669 ir = irow
670 END IF
671 120 CONTINUE
672 130 CONTINUE
673 140 CONTINUE
674*
675 jku = uub
676 DO 170 jkl = 1, llb
677*
678* Transform from bandwidth JKL-1, JKU to JKL, JKU
679*
680 DO 160 jc = 1, min( n+jkl, m ) + jku - 1
681 extra = zero
682 angle = twopi*slarnd( 1, iseed )
683 c = cos( angle )
684 s = sin( angle )
685 irow = max( 1, jc-jku )
686 IF( jc.LT.n ) THEN
687 il = min( m, jc+jkl ) + 1 - irow
688 CALL slarot( .false., jc.GT.jku, .false., il,
689 $ c,
690 $ s, a( irow-iskew*jc+ioffst, jc ),
691 $ ilda, extra, dummy )
692 END IF
693*
694* Chase "EXTRA" back up
695*
696 ic = jc
697 ir = irow
698 DO 150 jch = jc - jku, 1, -jkl - jku
699 IF( ic.LT.n ) THEN
700 CALL slartg( a( ir+1-iskew*( ic+1 )+ioffst,
701 $ ic+1 ), extra, c, s, dummy )
702 END IF
703 icol = max( 1, jch-jkl )
704 il = ic + 2 - icol
705 temp = zero
706 iltemp = jch.GT.jkl
707 CALL slarot( .true., iltemp, .true., il, c,
708 $ -s,
709 $ a( ir-iskew*icol+ioffst, icol ),
710 $ ilda, temp, extra )
711 IF( iltemp ) THEN
712 CALL slartg( a( ir+1-iskew*( icol+1 )+ioffst,
713 $ icol+1 ), temp, c, s, dummy )
714 irow = max( 1, jch-jkl-jku )
715 il = ir + 2 - irow
716 extra = zero
717 CALL slarot( .false., jch.GT.jkl+jku,
718 $ .true.,
719 $ il, c, -s, a( irow-iskew*icol+
720 $ ioffst, icol ), ilda, extra,
721 $ temp )
722 ic = icol
723 ir = irow
724 END IF
725 150 CONTINUE
726 160 CONTINUE
727 170 CONTINUE
728*
729 ELSE
730*
731* Bottom-Up -- Start at the bottom right.
732*
733 jkl = 0
734 DO 200 jku = 1, uub
735*
736* Transform from bandwidth JKL, JKU-1 to JKL, JKU
737*
738* First row actually rotated is M
739* First column actually rotated is MIN( M+JKU, N )
740*
741 iendch = min( m, n+jkl ) - 1
742 DO 190 jc = min( m+jku, n ) - 1, 1 - jkl, -1
743 extra = zero
744 angle = twopi*slarnd( 1, iseed )
745 c = cos( angle )
746 s = sin( angle )
747 irow = max( 1, jc-jku+1 )
748 IF( jc.GT.0 ) THEN
749 il = min( m, jc+jkl+1 ) + 1 - irow
750 CALL slarot( .false., .false., jc+jkl.LT.m,
751 $ il,
752 $ c, s, a( irow-iskew*jc+ioffst,
753 $ jc ), ilda, dummy, extra )
754 END IF
755*
756* Chase "EXTRA" back down
757*
758 ic = jc
759 DO 180 jch = jc + jkl, iendch, jkl + jku
760 ilextr = ic.GT.0
761 IF( ilextr ) THEN
762 CALL slartg( a( jch-iskew*ic+ioffst, ic ),
763 $ extra, c, s, dummy )
764 END IF
765 ic = max( 1, ic )
766 icol = min( n-1, jch+jku )
767 iltemp = jch + jku.LT.n
768 temp = zero
769 CALL slarot( .true., ilextr, iltemp,
770 $ icol+2-ic,
771 $ c, s, a( jch-iskew*ic+ioffst, ic ),
772 $ ilda, extra, temp )
773 IF( iltemp ) THEN
774 CALL slartg( a( jch-iskew*icol+ioffst,
775 $ icol ), temp, c, s, dummy )
776 il = min( iendch, jch+jkl+jku ) + 2 - jch
777 extra = zero
778 CALL slarot( .false., .true.,
779 $ jch+jkl+jku.LE.iendch, il, c, s,
780 $ a( jch-iskew*icol+ioffst,
781 $ icol ), ilda, temp, extra )
782 ic = icol
783 END IF
784 180 CONTINUE
785 190 CONTINUE
786 200 CONTINUE
787*
788 jku = uub
789 DO 230 jkl = 1, llb
790*
791* Transform from bandwidth JKL-1, JKU to JKL, JKU
792*
793* First row actually rotated is MIN( N+JKL, M )
794* First column actually rotated is N
795*
796 iendch = min( n, m+jku ) - 1
797 DO 220 jr = min( n+jkl, m ) - 1, 1 - jku, -1
798 extra = zero
799 angle = twopi*slarnd( 1, iseed )
800 c = cos( angle )
801 s = sin( angle )
802 icol = max( 1, jr-jkl+1 )
803 IF( jr.GT.0 ) THEN
804 il = min( n, jr+jku+1 ) + 1 - icol
805 CALL slarot( .true., .false., jr+jku.LT.n,
806 $ il,
807 $ c, s, a( jr-iskew*icol+ioffst,
808 $ icol ), ilda, dummy, extra )
809 END IF
810*
811* Chase "EXTRA" back down
812*
813 ir = jr
814 DO 210 jch = jr + jku, iendch, jkl + jku
815 ilextr = ir.GT.0
816 IF( ilextr ) THEN
817 CALL slartg( a( ir-iskew*jch+ioffst, jch ),
818 $ extra, c, s, dummy )
819 END IF
820 ir = max( 1, ir )
821 irow = min( m-1, jch+jkl )
822 iltemp = jch + jkl.LT.m
823 temp = zero
824 CALL slarot( .false., ilextr, iltemp,
825 $ irow+2-ir,
826 $ c, s, a( ir-iskew*jch+ioffst,
827 $ jch ), ilda, extra, temp )
828 IF( iltemp ) THEN
829 CALL slartg( a( irow-iskew*jch+ioffst, jch ),
830 $ temp, c, s, dummy )
831 il = min( iendch, jch+jkl+jku ) + 2 - jch
832 extra = zero
833 CALL slarot( .true., .true.,
834 $ jch+jkl+jku.LE.iendch, il, c, s,
835 $ a( irow-iskew*jch+ioffst, jch ),
836 $ ilda, temp, extra )
837 ir = irow
838 END IF
839 210 CONTINUE
840 220 CONTINUE
841 230 CONTINUE
842 END IF
843*
844 ELSE
845*
846* Symmetric -- A = U D U'
847*
848 ipackg = ipack
849 ioffg = ioffst
850*
851 IF( topdwn ) THEN
852*
853* Top-Down -- Generate Upper triangle only
854*
855 IF( ipack.GE.5 ) THEN
856 ipackg = 6
857 ioffg = uub + 1
858 ELSE
859 ipackg = 1
860 END IF
861 CALL scopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
862*
863 DO 260 k = 1, uub
864 DO 250 jc = 1, n - 1
865 irow = max( 1, jc-k )
866 il = min( jc+1, k+2 )
867 extra = zero
868 temp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
869 angle = twopi*slarnd( 1, iseed )
870 c = cos( angle )
871 s = sin( angle )
872 CALL slarot( .false., jc.GT.k, .true., il, c, s,
873 $ a( irow-iskew*jc+ioffg, jc ), ilda,
874 $ extra, temp )
875 CALL slarot( .true., .true., .false.,
876 $ min( k, n-jc )+1, c, s,
877 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
878 $ temp, dummy )
879*
880* Chase EXTRA back up the matrix
881*
882 icol = jc
883 DO 240 jch = jc - k, 1, -k
884 CALL slartg( a( jch+1-iskew*( icol+1 )+ioffg,
885 $ icol+1 ), extra, c, s, dummy )
886 temp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
887 CALL slarot( .true., .true., .true., k+2, c,
888 $ -s,
889 $ a( ( 1-iskew )*jch+ioffg, jch ),
890 $ ilda, temp, extra )
891 irow = max( 1, jch-k )
892 il = min( jch+1, k+2 )
893 extra = zero
894 CALL slarot( .false., jch.GT.k, .true., il,
895 $ c,
896 $ -s, a( irow-iskew*jch+ioffg, jch ),
897 $ ilda, extra, temp )
898 icol = jch
899 240 CONTINUE
900 250 CONTINUE
901 260 CONTINUE
902*
903* If we need lower triangle, copy from upper. Note that
904* the order of copying is chosen to work for 'q' -> 'b'
905*
906 IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
907 DO 280 jc = 1, n
908 irow = ioffst - iskew*jc
909 DO 270 jr = jc, min( n, jc+uub )
910 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
911 270 CONTINUE
912 280 CONTINUE
913 IF( ipack.EQ.5 ) THEN
914 DO 300 jc = n - uub + 1, n
915 DO 290 jr = n + 2 - jc, uub + 1
916 a( jr, jc ) = zero
917 290 CONTINUE
918 300 CONTINUE
919 END IF
920 IF( ipackg.EQ.6 ) THEN
921 ipackg = ipack
922 ELSE
923 ipackg = 0
924 END IF
925 END IF
926 ELSE
927*
928* Bottom-Up -- Generate Lower triangle only
929*
930 IF( ipack.GE.5 ) THEN
931 ipackg = 5
932 IF( ipack.EQ.6 )
933 $ ioffg = 1
934 ELSE
935 ipackg = 2
936 END IF
937 CALL scopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
938*
939 DO 330 k = 1, uub
940 DO 320 jc = n - 1, 1, -1
941 il = min( n+1-jc, k+2 )
942 extra = zero
943 temp = a( 1+( 1-iskew )*jc+ioffg, jc )
944 angle = twopi*slarnd( 1, iseed )
945 c = cos( angle )
946 s = -sin( angle )
947 CALL slarot( .false., .true., n-jc.GT.k, il, c,
948 $ s,
949 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
950 $ temp, extra )
951 icol = max( 1, jc-k+1 )
952 CALL slarot( .true., .false., .true., jc+2-icol,
953 $ c,
954 $ s, a( jc-iskew*icol+ioffg, icol ),
955 $ ilda, dummy, temp )
956*
957* Chase EXTRA back down the matrix
958*
959 icol = jc
960 DO 310 jch = jc + k, n - 1, k
961 CALL slartg( a( jch-iskew*icol+ioffg, icol ),
962 $ extra, c, s, dummy )
963 temp = a( 1+( 1-iskew )*jch+ioffg, jch )
964 CALL slarot( .true., .true., .true., k+2, c,
965 $ s,
966 $ a( jch-iskew*icol+ioffg, icol ),
967 $ ilda, extra, temp )
968 il = min( n+1-jch, k+2 )
969 extra = zero
970 CALL slarot( .false., .true., n-jch.GT.k, il,
971 $ c,
972 $ s, a( ( 1-iskew )*jch+ioffg, jch ),
973 $ ilda, temp, extra )
974 icol = jch
975 310 CONTINUE
976 320 CONTINUE
977 330 CONTINUE
978*
979* If we need upper triangle, copy from lower. Note that
980* the order of copying is chosen to work for 'b' -> 'q'
981*
982 IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
983 DO 350 jc = n, 1, -1
984 irow = ioffst - iskew*jc
985 DO 340 jr = jc, max( 1, jc-uub ), -1
986 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
987 340 CONTINUE
988 350 CONTINUE
989 IF( ipack.EQ.6 ) THEN
990 DO 370 jc = 1, uub
991 DO 360 jr = 1, uub + 1 - jc
992 a( jr, jc ) = zero
993 360 CONTINUE
994 370 CONTINUE
995 END IF
996 IF( ipackg.EQ.5 ) THEN
997 ipackg = ipack
998 ELSE
999 ipackg = 0
1000 END IF
1001 END IF
1002 END IF
1003 END IF
1004*
1005 ELSE
1006*
1007* 4) Generate Banded Matrix by first
1008* Rotating by random Unitary matrices,
1009* then reducing the bandwidth using Householder
1010* transformations.
1011*
1012* Note: we should get here only if LDA .ge. N
1013*
1014 IF( isym.EQ.1 ) THEN
1015*
1016* Non-symmetric -- A = U D V
1017*
1018 CALL slagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1019 $ iinfo )
1020 ELSE
1021*
1022* Symmetric -- A = U D U'
1023*
1024 CALL slagsy( m, llb, d, a, lda, iseed, work, iinfo )
1025*
1026 END IF
1027 IF( iinfo.NE.0 ) THEN
1028 info = 3
1029 RETURN
1030 END IF
1031 END IF
1032*
1033* 5) Pack the matrix
1034*
1035 IF( ipack.NE.ipackg ) THEN
1036 IF( ipack.EQ.1 ) THEN
1037*
1038* 'U' -- Upper triangular, not packed
1039*
1040 DO 390 j = 1, m
1041 DO 380 i = j + 1, m
1042 a( i, j ) = zero
1043 380 CONTINUE
1044 390 CONTINUE
1045*
1046 ELSE IF( ipack.EQ.2 ) THEN
1047*
1048* 'L' -- Lower triangular, not packed
1049*
1050 DO 410 j = 2, m
1051 DO 400 i = 1, j - 1
1052 a( i, j ) = zero
1053 400 CONTINUE
1054 410 CONTINUE
1055*
1056 ELSE IF( ipack.EQ.3 ) THEN
1057*
1058* 'C' -- Upper triangle packed Columnwise.
1059*
1060 icol = 1
1061 irow = 0
1062 DO 430 j = 1, m
1063 DO 420 i = 1, j
1064 irow = irow + 1
1065 IF( irow.GT.lda ) THEN
1066 irow = 1
1067 icol = icol + 1
1068 END IF
1069 a( irow, icol ) = a( i, j )
1070 420 CONTINUE
1071 430 CONTINUE
1072*
1073 ELSE IF( ipack.EQ.4 ) THEN
1074*
1075* 'R' -- Lower triangle packed Columnwise.
1076*
1077 icol = 1
1078 irow = 0
1079 DO 450 j = 1, m
1080 DO 440 i = j, m
1081 irow = irow + 1
1082 IF( irow.GT.lda ) THEN
1083 irow = 1
1084 icol = icol + 1
1085 END IF
1086 a( irow, icol ) = a( i, j )
1087 440 CONTINUE
1088 450 CONTINUE
1089*
1090 ELSE IF( ipack.GE.5 ) THEN
1091*
1092* 'B' -- The lower triangle is packed as a band matrix.
1093* 'Q' -- The upper triangle is packed as a band matrix.
1094* 'Z' -- The whole matrix is packed as a band matrix.
1095*
1096 IF( ipack.EQ.5 )
1097 $ uub = 0
1098 IF( ipack.EQ.6 )
1099 $ llb = 0
1100*
1101 DO 470 j = 1, uub
1102 DO 460 i = min( j+llb, m ), 1, -1
1103 a( i-j+uub+1, j ) = a( i, j )
1104 460 CONTINUE
1105 470 CONTINUE
1106*
1107 DO 490 j = uub + 2, n
1108 DO 480 i = j - uub, min( j+llb, m )
1109 a( i-j+uub+1, j ) = a( i, j )
1110 480 CONTINUE
1111 490 CONTINUE
1112 END IF
1113*
1114* If packed, zero out extraneous elements.
1115*
1116* Symmetric/Triangular Packed --
1117* zero out everything after A(IROW,ICOL)
1118*
1119 IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1120 DO 510 jc = icol, m
1121 DO 500 jr = irow + 1, lda
1122 a( jr, jc ) = zero
1123 500 CONTINUE
1124 irow = 0
1125 510 CONTINUE
1126*
1127 ELSE IF( ipack.GE.5 ) THEN
1128*
1129* Packed Band --
1130* 1st row is now in A( UUB+2-j, j), zero above it
1131* m-th row is now in A( M+UUB-j,j), zero below it
1132* last non-zero diagonal is now in A( UUB+LLB+1,j ),
1133* zero below it, too.
1134*
1135 ir1 = uub + llb + 2
1136 ir2 = uub + m + 2
1137 DO 540 jc = 1, n
1138 DO 520 jr = 1, uub + 1 - jc
1139 a( jr, jc ) = zero
1140 520 CONTINUE
1141 DO 530 jr = max( 1, min( ir1, ir2-jc ) ), lda
1142 a( jr, jc ) = zero
1143 530 CONTINUE
1144 540 CONTINUE
1145 END IF
1146 END IF
1147*
1148 RETURN
1149*
1150* End of SLATMT
1151*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine slagge(m, n, kl, ku, d, a, lda, iseed, work, info)
SLAGGE
Definition slagge.f:113
subroutine slagsy(n, k, d, a, lda, iseed, work, info)
SLAGSY
Definition slagsy.f:101
real function slarnd(idist, iseed)
SLARND
Definition slarnd.f:73
subroutine slarot(lrows, lleft, lright, nl, c, s, a, lda, xleft, xright)
SLAROT
Definition slarot.f:227
subroutine slatm7(mode, cond, irsign, idist, iseed, d, n, rank, info)
SLATM7
Definition slatm7.f:122
Here is the call graph for this function:
Here is the caller graph for this function: