LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ 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, "chasing"
              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, c,
635  $ s, a( jr-iskew*icol+ioffst, icol ),
636  $ ilda, extra, dummy )
637  END IF
638 *
639 * Chase "EXTRA" back up
640 *
641  ir = jr
642  ic = icol
643  DO 120 jch = jr - jkl, 1, -jkl - jku
644  IF( ir.LT.m ) THEN
645  CALL slartg( a( ir+1-iskew*( ic+1 )+ioffst,
646  $ ic+1 ), extra, c, s, dummy )
647  END IF
648  irow = max( 1, jch-jku )
649  il = ir + 2 - irow
650  temp = zero
651  iltemp = jch.GT.jku
652  CALL slarot( .false., iltemp, .true., il, c, -s,
653  $ a( irow-iskew*ic+ioffst, ic ),
654  $ ilda, temp, extra )
655  IF( iltemp ) THEN
656  CALL slartg( a( irow+1-iskew*( ic+1 )+ioffst,
657  $ ic+1 ), temp, c, s, dummy )
658  icol = max( 1, jch-jku-jkl )
659  il = ic + 2 - icol
660  extra = zero
661  CALL slarot( .true., jch.GT.jku+jkl, .true.,
662  $ il, c, -s, a( irow-iskew*icol+
663  $ ioffst, icol ), ilda, extra,
664  $ temp )
665  ic = icol
666  ir = irow
667  END IF
668  120 CONTINUE
669  130 CONTINUE
670  140 CONTINUE
671 *
672  jku = uub
673  DO 170 jkl = 1, llb
674 *
675 * Transform from bandwidth JKL-1, JKU to JKL, JKU
676 *
677  DO 160 jc = 1, min( n+jkl, m ) + jku - 1
678  extra = zero
679  angle = twopi*slarnd( 1, iseed )
680  c = cos( angle )
681  s = sin( angle )
682  irow = max( 1, jc-jku )
683  IF( jc.LT.n ) THEN
684  il = min( m, jc+jkl ) + 1 - irow
685  CALL slarot( .false., jc.GT.jku, .false., il, c,
686  $ s, a( irow-iskew*jc+ioffst, jc ),
687  $ ilda, extra, dummy )
688  END IF
689 *
690 * Chase "EXTRA" back up
691 *
692  ic = jc
693  ir = irow
694  DO 150 jch = jc - jku, 1, -jkl - jku
695  IF( ic.LT.n ) THEN
696  CALL slartg( a( ir+1-iskew*( ic+1 )+ioffst,
697  $ ic+1 ), extra, c, s, dummy )
698  END IF
699  icol = max( 1, jch-jkl )
700  il = ic + 2 - icol
701  temp = zero
702  iltemp = jch.GT.jkl
703  CALL slarot( .true., iltemp, .true., il, c, -s,
704  $ a( ir-iskew*icol+ioffst, icol ),
705  $ ilda, temp, extra )
706  IF( iltemp ) THEN
707  CALL slartg( a( ir+1-iskew*( icol+1 )+ioffst,
708  $ icol+1 ), temp, c, s, dummy )
709  irow = max( 1, jch-jkl-jku )
710  il = ir + 2 - irow
711  extra = zero
712  CALL slarot( .false., jch.GT.jkl+jku, .true.,
713  $ il, c, -s, a( irow-iskew*icol+
714  $ ioffst, icol ), ilda, extra,
715  $ temp )
716  ic = icol
717  ir = irow
718  END IF
719  150 CONTINUE
720  160 CONTINUE
721  170 CONTINUE
722 *
723  ELSE
724 *
725 * Bottom-Up -- Start at the bottom right.
726 *
727  jkl = 0
728  DO 200 jku = 1, uub
729 *
730 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
731 *
732 * First row actually rotated is M
733 * First column actually rotated is MIN( M+JKU, N )
734 *
735  iendch = min( m, n+jkl ) - 1
736  DO 190 jc = min( m+jku, n ) - 1, 1 - jkl, -1
737  extra = zero
738  angle = twopi*slarnd( 1, iseed )
739  c = cos( angle )
740  s = sin( angle )
741  irow = max( 1, jc-jku+1 )
742  IF( jc.GT.0 ) THEN
743  il = min( m, jc+jkl+1 ) + 1 - irow
744  CALL slarot( .false., .false., jc+jkl.LT.m, il,
745  $ c, s, a( irow-iskew*jc+ioffst,
746  $ jc ), ilda, dummy, extra )
747  END IF
748 *
749 * Chase "EXTRA" back down
750 *
751  ic = jc
752  DO 180 jch = jc + jkl, iendch, jkl + jku
753  ilextr = ic.GT.0
754  IF( ilextr ) THEN
755  CALL slartg( a( jch-iskew*ic+ioffst, ic ),
756  $ extra, c, s, dummy )
757  END IF
758  ic = max( 1, ic )
759  icol = min( n-1, jch+jku )
760  iltemp = jch + jku.LT.n
761  temp = zero
762  CALL slarot( .true., ilextr, iltemp, icol+2-ic,
763  $ c, s, a( jch-iskew*ic+ioffst, ic ),
764  $ ilda, extra, temp )
765  IF( iltemp ) THEN
766  CALL slartg( a( jch-iskew*icol+ioffst,
767  $ icol ), temp, c, s, dummy )
768  il = min( iendch, jch+jkl+jku ) + 2 - jch
769  extra = zero
770  CALL slarot( .false., .true.,
771  $ jch+jkl+jku.LE.iendch, il, c, s,
772  $ a( jch-iskew*icol+ioffst,
773  $ icol ), ilda, temp, extra )
774  ic = icol
775  END IF
776  180 CONTINUE
777  190 CONTINUE
778  200 CONTINUE
779 *
780  jku = uub
781  DO 230 jkl = 1, llb
782 *
783 * Transform from bandwidth JKL-1, JKU to JKL, JKU
784 *
785 * First row actually rotated is MIN( N+JKL, M )
786 * First column actually rotated is N
787 *
788  iendch = min( n, m+jku ) - 1
789  DO 220 jr = min( n+jkl, m ) - 1, 1 - jku, -1
790  extra = zero
791  angle = twopi*slarnd( 1, iseed )
792  c = cos( angle )
793  s = sin( angle )
794  icol = max( 1, jr-jkl+1 )
795  IF( jr.GT.0 ) THEN
796  il = min( n, jr+jku+1 ) + 1 - icol
797  CALL slarot( .true., .false., jr+jku.LT.n, il,
798  $ c, s, a( jr-iskew*icol+ioffst,
799  $ icol ), ilda, dummy, extra )
800  END IF
801 *
802 * Chase "EXTRA" back down
803 *
804  ir = jr
805  DO 210 jch = jr + jku, iendch, jkl + jku
806  ilextr = ir.GT.0
807  IF( ilextr ) THEN
808  CALL slartg( a( ir-iskew*jch+ioffst, jch ),
809  $ extra, c, s, dummy )
810  END IF
811  ir = max( 1, ir )
812  irow = min( m-1, jch+jkl )
813  iltemp = jch + jkl.LT.m
814  temp = zero
815  CALL slarot( .false., ilextr, iltemp, irow+2-ir,
816  $ c, s, a( ir-iskew*jch+ioffst,
817  $ jch ), ilda, extra, temp )
818  IF( iltemp ) THEN
819  CALL slartg( a( irow-iskew*jch+ioffst, jch ),
820  $ temp, c, s, dummy )
821  il = min( iendch, jch+jkl+jku ) + 2 - jch
822  extra = zero
823  CALL slarot( .true., .true.,
824  $ jch+jkl+jku.LE.iendch, il, c, s,
825  $ a( irow-iskew*jch+ioffst, jch ),
826  $ ilda, temp, extra )
827  ir = irow
828  END IF
829  210 CONTINUE
830  220 CONTINUE
831  230 CONTINUE
832  END IF
833 *
834  ELSE
835 *
836 * Symmetric -- A = U D U'
837 *
838  ipackg = ipack
839  ioffg = ioffst
840 *
841  IF( topdwn ) THEN
842 *
843 * Top-Down -- Generate Upper triangle only
844 *
845  IF( ipack.GE.5 ) THEN
846  ipackg = 6
847  ioffg = uub + 1
848  ELSE
849  ipackg = 1
850  END IF
851  CALL scopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
852 *
853  DO 260 k = 1, uub
854  DO 250 jc = 1, n - 1
855  irow = max( 1, jc-k )
856  il = min( jc+1, k+2 )
857  extra = zero
858  temp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
859  angle = twopi*slarnd( 1, iseed )
860  c = cos( angle )
861  s = sin( angle )
862  CALL slarot( .false., jc.GT.k, .true., il, c, s,
863  $ a( irow-iskew*jc+ioffg, jc ), ilda,
864  $ extra, temp )
865  CALL slarot( .true., .true., .false.,
866  $ min( k, n-jc )+1, c, s,
867  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
868  $ temp, dummy )
869 *
870 * Chase EXTRA back up the matrix
871 *
872  icol = jc
873  DO 240 jch = jc - k, 1, -k
874  CALL slartg( a( jch+1-iskew*( icol+1 )+ioffg,
875  $ icol+1 ), extra, c, s, dummy )
876  temp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
877  CALL slarot( .true., .true., .true., k+2, c, -s,
878  $ a( ( 1-iskew )*jch+ioffg, jch ),
879  $ ilda, temp, extra )
880  irow = max( 1, jch-k )
881  il = min( jch+1, k+2 )
882  extra = zero
883  CALL slarot( .false., jch.GT.k, .true., il, c,
884  $ -s, a( irow-iskew*jch+ioffg, jch ),
885  $ ilda, extra, temp )
886  icol = jch
887  240 CONTINUE
888  250 CONTINUE
889  260 CONTINUE
890 *
891 * If we need lower triangle, copy from upper. Note that
892 * the order of copying is chosen to work for 'q' -> 'b'
893 *
894  IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
895  DO 280 jc = 1, n
896  irow = ioffst - iskew*jc
897  DO 270 jr = jc, min( n, jc+uub )
898  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
899  270 CONTINUE
900  280 CONTINUE
901  IF( ipack.EQ.5 ) THEN
902  DO 300 jc = n - uub + 1, n
903  DO 290 jr = n + 2 - jc, uub + 1
904  a( jr, jc ) = zero
905  290 CONTINUE
906  300 CONTINUE
907  END IF
908  IF( ipackg.EQ.6 ) THEN
909  ipackg = ipack
910  ELSE
911  ipackg = 0
912  END IF
913  END IF
914  ELSE
915 *
916 * Bottom-Up -- Generate Lower triangle only
917 *
918  IF( ipack.GE.5 ) THEN
919  ipackg = 5
920  IF( ipack.EQ.6 )
921  $ ioffg = 1
922  ELSE
923  ipackg = 2
924  END IF
925  CALL scopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
926 *
927  DO 330 k = 1, uub
928  DO 320 jc = n - 1, 1, -1
929  il = min( n+1-jc, k+2 )
930  extra = zero
931  temp = a( 1+( 1-iskew )*jc+ioffg, jc )
932  angle = twopi*slarnd( 1, iseed )
933  c = cos( angle )
934  s = -sin( angle )
935  CALL slarot( .false., .true., n-jc.GT.k, il, c, s,
936  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
937  $ temp, extra )
938  icol = max( 1, jc-k+1 )
939  CALL slarot( .true., .false., .true., jc+2-icol, c,
940  $ s, a( jc-iskew*icol+ioffg, icol ),
941  $ ilda, dummy, temp )
942 *
943 * Chase EXTRA back down the matrix
944 *
945  icol = jc
946  DO 310 jch = jc + k, n - 1, k
947  CALL slartg( a( jch-iskew*icol+ioffg, icol ),
948  $ extra, c, s, dummy )
949  temp = a( 1+( 1-iskew )*jch+ioffg, jch )
950  CALL slarot( .true., .true., .true., k+2, c, s,
951  $ a( jch-iskew*icol+ioffg, icol ),
952  $ ilda, extra, temp )
953  il = min( n+1-jch, k+2 )
954  extra = zero
955  CALL slarot( .false., .true., n-jch.GT.k, il, c,
956  $ s, a( ( 1-iskew )*jch+ioffg, jch ),
957  $ ilda, temp, extra )
958  icol = jch
959  310 CONTINUE
960  320 CONTINUE
961  330 CONTINUE
962 *
963 * If we need upper triangle, copy from lower. Note that
964 * the order of copying is chosen to work for 'b' -> 'q'
965 *
966  IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
967  DO 350 jc = n, 1, -1
968  irow = ioffst - iskew*jc
969  DO 340 jr = jc, max( 1, jc-uub ), -1
970  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
971  340 CONTINUE
972  350 CONTINUE
973  IF( ipack.EQ.6 ) THEN
974  DO 370 jc = 1, uub
975  DO 360 jr = 1, uub + 1 - jc
976  a( jr, jc ) = zero
977  360 CONTINUE
978  370 CONTINUE
979  END IF
980  IF( ipackg.EQ.5 ) THEN
981  ipackg = ipack
982  ELSE
983  ipackg = 0
984  END IF
985  END IF
986  END IF
987  END IF
988 *
989  ELSE
990 *
991 * 4) Generate Banded Matrix by first
992 * Rotating by random Unitary matrices,
993 * then reducing the bandwidth using Householder
994 * transformations.
995 *
996 * Note: we should get here only if LDA .ge. N
997 *
998  IF( isym.EQ.1 ) THEN
999 *
1000 * Non-symmetric -- A = U D V
1001 *
1002  CALL slagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1003  $ iinfo )
1004  ELSE
1005 *
1006 * Symmetric -- A = U D U'
1007 *
1008  CALL slagsy( m, llb, d, a, lda, iseed, work, iinfo )
1009 *
1010  END IF
1011  IF( iinfo.NE.0 ) THEN
1012  info = 3
1013  RETURN
1014  END IF
1015  END IF
1016 *
1017 * 5) Pack the matrix
1018 *
1019  IF( ipack.NE.ipackg ) THEN
1020  IF( ipack.EQ.1 ) THEN
1021 *
1022 * 'U' -- Upper triangular, not packed
1023 *
1024  DO 390 j = 1, m
1025  DO 380 i = j + 1, m
1026  a( i, j ) = zero
1027  380 CONTINUE
1028  390 CONTINUE
1029 *
1030  ELSE IF( ipack.EQ.2 ) THEN
1031 *
1032 * 'L' -- Lower triangular, not packed
1033 *
1034  DO 410 j = 2, m
1035  DO 400 i = 1, j - 1
1036  a( i, j ) = zero
1037  400 CONTINUE
1038  410 CONTINUE
1039 *
1040  ELSE IF( ipack.EQ.3 ) THEN
1041 *
1042 * 'C' -- Upper triangle packed Columnwise.
1043 *
1044  icol = 1
1045  irow = 0
1046  DO 430 j = 1, m
1047  DO 420 i = 1, j
1048  irow = irow + 1
1049  IF( irow.GT.lda ) THEN
1050  irow = 1
1051  icol = icol + 1
1052  END IF
1053  a( irow, icol ) = a( i, j )
1054  420 CONTINUE
1055  430 CONTINUE
1056 *
1057  ELSE IF( ipack.EQ.4 ) THEN
1058 *
1059 * 'R' -- Lower triangle packed Columnwise.
1060 *
1061  icol = 1
1062  irow = 0
1063  DO 450 j = 1, m
1064  DO 440 i = j, m
1065  irow = irow + 1
1066  IF( irow.GT.lda ) THEN
1067  irow = 1
1068  icol = icol + 1
1069  END IF
1070  a( irow, icol ) = a( i, j )
1071  440 CONTINUE
1072  450 CONTINUE
1073 *
1074  ELSE IF( ipack.GE.5 ) THEN
1075 *
1076 * 'B' -- The lower triangle is packed as a band matrix.
1077 * 'Q' -- The upper triangle is packed as a band matrix.
1078 * 'Z' -- The whole matrix is packed as a band matrix.
1079 *
1080  IF( ipack.EQ.5 )
1081  $ uub = 0
1082  IF( ipack.EQ.6 )
1083  $ llb = 0
1084 *
1085  DO 470 j = 1, uub
1086  DO 460 i = min( j+llb, m ), 1, -1
1087  a( i-j+uub+1, j ) = a( i, j )
1088  460 CONTINUE
1089  470 CONTINUE
1090 *
1091  DO 490 j = uub + 2, n
1092  DO 480 i = j - uub, min( j+llb, m )
1093  a( i-j+uub+1, j ) = a( i, j )
1094  480 CONTINUE
1095  490 CONTINUE
1096  END IF
1097 *
1098 * If packed, zero out extraneous elements.
1099 *
1100 * Symmetric/Triangular Packed --
1101 * zero out everything after A(IROW,ICOL)
1102 *
1103  IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1104  DO 510 jc = icol, m
1105  DO 500 jr = irow + 1, lda
1106  a( jr, jc ) = zero
1107  500 CONTINUE
1108  irow = 0
1109  510 CONTINUE
1110 *
1111  ELSE IF( ipack.GE.5 ) THEN
1112 *
1113 * Packed Band --
1114 * 1st row is now in A( UUB+2-j, j), zero above it
1115 * m-th row is now in A( M+UUB-j,j), zero below it
1116 * last non-zero diagonal is now in A( UUB+LLB+1,j ),
1117 * zero below it, too.
1118 *
1119  ir1 = uub + llb + 2
1120  ir2 = uub + m + 2
1121  DO 540 jc = 1, n
1122  DO 520 jr = 1, uub + 1 - jc
1123  a( jr, jc ) = zero
1124  520 CONTINUE
1125  DO 530 jr = max( 1, min( ir1, ir2-jc ) ), lda
1126  a( jr, jc ) = zero
1127  530 CONTINUE
1128  540 CONTINUE
1129  END IF
1130  END IF
1131 *
1132  RETURN
1133 *
1134 * End of SLATMT
1135 *
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:110
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:97
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine slatm7(MODE, COND, IRSIGN, IDIST, ISEED, D, N, RANK, INFO)
SLATM7
Definition: slatm7.f:122
subroutine slarot(LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
SLAROT
Definition: slarot.f:226
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 slagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
SLAGGE
Definition: slagge.f:113
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: