LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlatms()

subroutine dlatms ( integer  M,
integer  N,
character  DIST,
integer, dimension( 4 )  ISEED,
character  SYM,
double precision, dimension( * )  D,
integer  MODE,
double precision  COND,
double precision  DMAX,
integer  KL,
integer  KU,
character  PACK,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLATMS

Purpose:
    DLATMS generates random matrices with specified singular values
    (or symmetric/hermitian with specified eigenvalues)
    for testing LAPACK programs.

    DLATMS 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 DLATMS
           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 DOUBLE PRECISION 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:N)=1.0/COND
           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
           MODE = 3 sets D(I)=COND**(-(I-1)/(N-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 DOUBLE PRECISION
           On entry, this is used as described under MODE above.
           If used, it must be >= 1. Not modified.
[in]DMAX
          DMAX is DOUBLE PRECISION
           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]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 DLATMS differ only in the PACK parameter,
           they will generate mathematically equivalent matrices.
           Not modified.
[in,out]A
          A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLATM1
            2  => Cannot scale to DMAX (max. sing. value is 0)
            3  => Error return from DLAGGE or SLAGSY
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 319 of file dlatms.f.

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