LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 323 of file dlatms.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: