zlatms.f

Go to the documentation of this file.

      SUBROUTINE zlatms( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
     $                   KL, KU, PACK, A, LDA, WORK, INFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          DIST, PACK, SYM
      INTEGER            INFO, KL, KU, LDA, M, MODE, N
      DOUBLE PRECISION   COND, DMAX
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 )
      DOUBLE PRECISION   D( * )
      COMPLEX*16         A( LDA, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*     ZLATMS generates random matrices with specified singular values
*     (or hermitian with specified eigenvalues)
*     for testing LAPACK programs.
*
*     ZLATMS 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 hermitian or 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 hermitian or 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 hermitian)
*          zero out lower half (if hermitian)
*          store the upper half columnwise (if hermitian or upper
*                triangular)
*          store the lower half columnwise (if hermitian or lower
*                triangular)
*          store the lower triangle in banded format (if hermitian or
*                lower triangular)
*          store the upper triangle in banded format (if hermitian 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.
*
*  Arguments
*  =========
*
*  M      - INTEGER
*           The number of rows of A. Not modified.
*
*  N      - INTEGER
*           The number of columns of A. N must equal M if the matrix
*           is symmetric or hermitian (i.e., if SYM is not 'N')
*           Not modified.
*
*  DIST   - 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.
*
*  ISEED  - 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 ZLATMS
*           to continue the same random number sequence.
*           Changed on exit.
*
*  SYM    - CHARACTER*1
*           If SYM='H', the generated matrix is hermitian, with
*             eigenvalues specified by D, COND, MODE, and DMAX; they
*             may be positive, negative, or zero.
*           If SYM='P', the generated matrix is hermitian, 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.
*           If SYM='S', the generated matrix is (complex) symmetric,
*             with singular values specified by D, COND, MODE, and
*             DMAX; they will not be negative.
*           Not modified.
*
*  D      - 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.
*
*  MODE   - 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='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.
*
*  COND   - DOUBLE PRECISION
*           On entry, this is used as described under MODE above.
*           If used, it must be >= 1. Not modified.
*
*  DMAX   - 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.
*
*  KL     - 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 or hermitian.
*           Not modified.
*
*  KU     - 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 or hermitian.
*           Not modified.
*
*  PACK   - CHARACTER*1
*           This specifies packing of matrix as follows:
*           'N' => no packing
*           'U' => zero out all subdiagonal entries (if symmetric
*                  or hermitian)
*           'L' => zero out all superdiagonal entries (if symmetric
*                  or hermitian)
*           'C' => store the upper triangle columnwise (only if the
*                  matrix is symmetric, hermitian, or upper triangular)
*           'R' => store the lower triangle columnwise (only if the
*                  matrix is symmetric, hermitian, or lower triangular)
*           'B' => store the lower triangle in band storage scheme
*                  (only if the matrix is symmetric, hermitian, or
*                  lower triangular)
*           'Q' => store the upper triangle in band storage scheme
*                  (only if the matrix is symmetric, hermitian, 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, HB, or TB     - use 'B' or 'Q'
*           PP, SP, HB, or TP     - use 'C' or 'R'
*
*           If two calls to ZLATMS differ only in the PACK parameter,
*           they will generate mathematically equivalent matrices.
*           Not modified.
*
*  A      - COMPLEX*16 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.
*
*  LDA    - 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.
*
*  WORK   - COMPLEX*16 array, dimension ( 3*MAX( N, M ) )
*           Workspace.
*           Modified.
*
*  INFO   - 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 is not 'N' and KU is 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 ZLAGGE, CLAGHE or CLAGSY
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d+0 )
      DOUBLE PRECISION   ONE
      parameter( one = 1.0d+0 )
      COMPLEX*16         CZERO
      parameter( czero = ( 0.0d+0, 0.0d+0 ) )
      DOUBLE PRECISION   TWOPI
      parameter( twopi = 6.2831853071795864769252867663d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            GIVENS, ILEXTR, ILTEMP, TOPDWN, ZSYM
      INTEGER            I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
     $                   ioffg, ioffst, ipack, ipackg, ir, ir1, ir2,
     $                   irow, irsign, iskew, isym, isympk, j, jc, jch,
     $                   jkl, jku, jr, k, llb, minlda, mnmin, mr, nc,
     $                   uub
      DOUBLE PRECISION   ALPHA, ANGLE, REALC, TEMP
      COMPLEX*16         C, CT, CTEMP, DUMMY, EXTRA, S, ST
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLARND
      COMPLEX*16         ZLARND
      EXTERNAL           lsame, dlarnd, zlarnd
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlatm1, dscal, xerbla, zlagge, zlaghe, zlagsy,
     $                   zlarot, zlartg, zlaset
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, cos, dble, dcmplx, dconjg, max, min, mod,
     $                   sin
*     ..
*     .. Executable Statements ..
*
*     1)      Decode and Test the input parameters.
*             Initialize flags & seed.
*
      info = 0
*
*     Quick return if possible
*
      IF( m.EQ.0 .OR. n.EQ.0 )
     $   RETURN
*
*     Decode DIST
*
      IF( lsame( dist, 'U' ) ) THEN
         idist = 1
      ELSE IF( lsame( dist, 'S' ) ) THEN
         idist = 2
      ELSE IF( lsame( dist, 'N' ) ) THEN
         idist = 3
      ELSE
         idist = -1
      END IF
*
*     Decode SYM
*
      IF( lsame( sym, 'N' ) ) THEN
         isym = 1
         irsign = 0
         zsym = .false.
      ELSE IF( lsame( sym, 'P' ) ) THEN
         isym = 2
         irsign = 0
         zsym = .false.
      ELSE IF( lsame( sym, 'S' ) ) THEN
         isym = 2
         irsign = 0
         zsym = .true.
      ELSE IF( lsame( sym, 'H' ) ) THEN
         isym = 2
         irsign = 1
         zsym = .false.
      ELSE
         isym = -1
      END IF
*
*     Decode PACK
*
      isympk = 0
      IF( lsame( pack, 'N' ) ) THEN
         ipack = 0
      ELSE IF( lsame( pack, 'U' ) ) THEN
         ipack = 1
         isympk = 1
      ELSE IF( lsame( pack, 'L' ) ) THEN
         ipack = 2
         isympk = 1
      ELSE IF( lsame( pack, 'C' ) ) THEN
         ipack = 3
         isympk = 2
      ELSE IF( lsame( pack, 'R' ) ) THEN
         ipack = 4
         isympk = 3
      ELSE IF( lsame( pack, 'B' ) ) THEN
         ipack = 5
         isympk = 3
      ELSE IF( lsame( pack, 'Q' ) ) THEN
         ipack = 6
         isympk = 2
      ELSE IF( lsame( pack, 'Z' ) ) THEN
         ipack = 7
      ELSE
         ipack = -1
      END IF
*
*     Set certain internal parameters
*
      mnmin = min( m, n )
      llb = min( kl, m-1 )
      uub = min( ku, n-1 )
      mr = min( m, n+llb )
      nc = min( n, m+uub )
      irow = 1
      icol = 1
      zsym = .false.
*
      IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
         minlda = uub + 1
      ELSE IF( ipack.EQ.7 ) THEN
         minlda = llb + uub + 1
      ELSE
         minlda = m
      END IF
*
*     Use Givens rotation method if bandwidth small enough,
*     or if LDA is too small to store the matrix unpacked.
*
      givens = .false.
      IF( isym.EQ.1 ) THEN
         IF( dble( llb+uub ).LT.0.3d0*dble( max( 1, mr+nc ) ) )
     $      givens = .true.
      ELSE
         IF( 2*llb.LT.m )
     $      givens = .true.
      END IF
      IF( lda.LT.m .AND. lda.GE.minlda )
     $   givens = .true.
*
*     Set INFO if an error
*
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( idist.EQ.-1 ) THEN
         info = -3
      ELSE IF( isym.EQ.-1 ) THEN
         info = -5
      ELSE IF( abs( mode ).GT.6 ) THEN
         info = -7
      ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
     $          THEN
         info = -8
      ELSE IF( kl.LT.0 ) THEN
         info = -10
      ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
         info = -11
      ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
     $         ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
     $         ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
     $         ( isympk.NE.0 .AND. m.NE.n ) ) THEN
         info = -12
      ELSE IF( lda.LT.max( 1, minlda ) ) THEN
         info = -14
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZLATMS', -info )
         RETURN
      END IF
*
*     Initialize random number generator
*
      DO 10 i = 1, 4
         iseed( i ) = mod( abs( iseed( i ) ), 4096 )
   10 CONTINUE
*
      IF( mod( iseed( 4 ), 2 ).NE.1 )
     $   iseed( 4 ) = iseed( 4 ) + 1
*
*     2)      Set up D  if indicated.
*
*             Compute D according to COND and MODE
*
      CALL dlatm1( mode, cond, irsign, idist, iseed, d, mnmin, iinfo )
      IF( iinfo.NE.0 ) THEN
         info = 1
         RETURN
      END IF
*
*     Choose Top-Down if D is (apparently) increasing,
*     Bottom-Up if D is (apparently) decreasing.
*
      IF( abs( d( 1 ) ).LE.abs( d( mnmin ) ) ) THEN
         topdwn = .true.
      ELSE
         topdwn = .false.
      END IF
*
      IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
*
*        Scale by DMAX
*
         temp = abs( d( 1 ) )
         DO 20 i = 2, mnmin
            temp = max( temp, abs( d( i ) ) )
   20    CONTINUE
*
         IF( temp.GT.zero ) THEN
            alpha = dmax / temp
         ELSE
            info = 2
            RETURN
         END IF
*
         CALL dscal( mnmin, alpha, d, 1 )
*
      END IF
*
      CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
*
*     3)      Generate Banded Matrix using Givens rotations.
*             Also the special case of UUB=LLB=0
*
*               Compute Addressing constants to cover all
*               storage formats.  Whether GE, HE, SY, GB, HB, or SB,
*               upper or lower triangle or both,
*               the (i,j)-th element is in
*               A( i - ISKEW*j + IOFFST, j )
*
      IF( ipack.GT.4 ) THEN
         ilda = lda - 1
         iskew = 1
         IF( ipack.GT.5 ) THEN
            ioffst = uub + 1
         ELSE
            ioffst = 1
         END IF
      ELSE
         ilda = lda
         iskew = 0
         ioffst = 0
      END IF
*
*     IPACKG is the format that the matrix is generated in. If this is
*     different from IPACK, then the matrix must be repacked at the
*     end.  It also signals how to compute the norm, for scaling.
*
      ipackg = 0
*
*     Diagonal Matrix -- We are done, unless it
*     is to be stored HP/SP/PP/TP (PACK='R' or 'C')
*
      IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
         DO 30 j = 1, mnmin
            a( ( 1-iskew )*j+ioffst, j ) = dcmplx( d( j ) )
   30    CONTINUE
*
         IF( ipack.LE.2 .OR. ipack.GE.5 )
     $      ipackg = ipack
*
      ELSE IF( givens ) THEN
*
*        Check whether to use Givens rotations,
*        Householder transformations, or nothing.
*
         IF( isym.EQ.1 ) THEN
*
*           Non-symmetric -- A = U D V
*
            IF( ipack.GT.4 ) THEN
               ipackg = ipack
            ELSE
               ipackg = 0
            END IF
*
            DO 40 j = 1, mnmin
               a( ( 1-iskew )*j+ioffst, j ) = dcmplx( d( j ) )
   40       CONTINUE
*
            IF( topdwn ) THEN
               jkl = 0
               DO 70 jku = 1, uub
*
*                 Transform from bandwidth JKL, JKU-1 to JKL, JKU
*
*                 Last row actually rotated is M
*                 Last column actually rotated is MIN( M+JKU, N )
*
                  DO 60 jr = 1, min( m+jku, n ) + jkl - 1
                     extra = czero
                     angle = twopi*dlarnd( 1, iseed )
                     c = cos( angle )*zlarnd( 5, iseed )
                     s = sin( angle )*zlarnd( 5, iseed )
                     icol = max( 1, jr-jkl )
                     IF( jr.LT.m ) THEN
                        il = min( n, jr+jku ) + 1 - icol
                        CALL zlarot( .true., jr.GT.jkl, .false., il, c,
     $                               s, a( jr-iskew*icol+ioffst, icol ),
     $                               ilda, extra, dummy )
                     END IF
*
*                    Chase "EXTRA" back up
*
                     ir = jr
                     ic = icol
                     DO 50 jch = jr - jkl, 1, -jkl - jku
                        IF( ir.LT.m ) THEN
                           CALL zlartg( a( ir+1-iskew*( ic+1 )+ioffst,
     $                                  ic+1 ), extra, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = dconjg( realc*dummy )
                           s = dconjg( -s*dummy )
                        END IF
                        irow = max( 1, jch-jku )
                        il = ir + 2 - irow
                        ctemp = czero
                        iltemp = jch.GT.jku
                        CALL zlarot( .false., iltemp, .true., il, c, s,
     $                               a( irow-iskew*ic+ioffst, ic ),
     $                               ilda, ctemp, extra )
                        IF( iltemp ) THEN
                           CALL zlartg( a( irow+1-iskew*( ic+1 )+ioffst,
     $                                  ic+1 ), ctemp, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = dconjg( realc*dummy )
                           s = dconjg( -s*dummy )
*
                           icol = max( 1, jch-jku-jkl )
                           il = ic + 2 - icol
                           extra = czero
                           CALL zlarot( .true., jch.GT.jku+jkl, .true.,
     $                                  il, c, s, a( irow-iskew*icol+
     $                                  ioffst, icol ), ilda, extra,
     $                                  ctemp )
                           ic = icol
                           ir = irow
                        END IF
   50                CONTINUE
   60             CONTINUE
   70          CONTINUE
*
               jku = uub
               DO 100 jkl = 1, llb
*
*                 Transform from bandwidth JKL-1, JKU to JKL, JKU
*
                  DO 90 jc = 1, min( n+jkl, m ) + jku - 1
                     extra = czero
                     angle = twopi*dlarnd( 1, iseed )
                     c = cos( angle )*zlarnd( 5, iseed )
                     s = sin( angle )*zlarnd( 5, iseed )
                     irow = max( 1, jc-jku )
                     IF( jc.LT.n ) THEN
                        il = min( m, jc+jkl ) + 1 - irow
                        CALL zlarot( .false., jc.GT.jku, .false., il, c,
     $                               s, a( irow-iskew*jc+ioffst, jc ),
     $                               ilda, extra, dummy )
                     END IF
*
*                    Chase "EXTRA" back up
*
                     ic = jc
                     ir = irow
                     DO 80 jch = jc - jku, 1, -jkl - jku
                        IF( ic.LT.n ) THEN
                           CALL zlartg( a( ir+1-iskew*( ic+1 )+ioffst,
     $                                  ic+1 ), extra, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = dconjg( realc*dummy )
                           s = dconjg( -s*dummy )
                        END IF
                        icol = max( 1, jch-jkl )
                        il = ic + 2 - icol
                        ctemp = czero
                        iltemp = jch.GT.jkl
                        CALL zlarot( .true., iltemp, .true., il, c, s,
     $                               a( ir-iskew*icol+ioffst, icol ),
     $                               ilda, ctemp, extra )
                        IF( iltemp ) THEN
                           CALL zlartg( a( ir+1-iskew*( icol+1 )+ioffst,
     $                                  icol+1 ), ctemp, realc, s,
     $                                  dummy )
                           dummy = zlarnd( 5, iseed )
                           c = dconjg( realc*dummy )
                           s = dconjg( -s*dummy )
                           irow = max( 1, jch-jkl-jku )
                           il = ir + 2 - irow
                           extra = czero
                           CALL zlarot( .false., jch.GT.jkl+jku, .true.,
     $                                  il, c, s, a( irow-iskew*icol+
     $                                  ioffst, icol ), ilda, extra,
     $                                  ctemp )
                           ic = icol
                           ir = irow
                        END IF
   80                CONTINUE
   90             CONTINUE
  100          CONTINUE
*
            ELSE
*
*              Bottom-Up -- Start at the bottom right.
*
               jkl = 0
               DO 130 jku = 1, uub
*
*                 Transform from bandwidth JKL, JKU-1 to JKL, JKU
*
*                 First row actually rotated is M
*                 First column actually rotated is MIN( M+JKU, N )
*
                  iendch = min( m, n+jkl ) - 1
                  DO 120 jc = min( m+jku, n ) - 1, 1 - jkl, -1
                     extra = czero
                     angle = twopi*dlarnd( 1, iseed )
                     c = cos( angle )*zlarnd( 5, iseed )
                     s = sin( angle )*zlarnd( 5, iseed )
                     irow = max( 1, jc-jku+1 )
                     IF( jc.GT.0 ) THEN
                        il = min( m, jc+jkl+1 ) + 1 - irow
                        CALL zlarot( .false., .false., jc+jkl.LT.m, il,
     $                               c, s, a( irow-iskew*jc+ioffst,
     $                               jc ), ilda, dummy, extra )
                     END IF
*
*                    Chase "EXTRA" back down
*
                     ic = jc
                     DO 110 jch = jc + jkl, iendch, jkl + jku
                        ilextr = ic.GT.0
                        IF( ilextr ) THEN
                           CALL zlartg( a( jch-iskew*ic+ioffst, ic ),
     $                                  extra, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = realc*dummy
                           s = s*dummy
                        END IF
                        ic = max( 1, ic )
                        icol = min( n-1, jch+jku )
                        iltemp = jch + jku.LT.n
                        ctemp = czero
                        CALL zlarot( .true., ilextr, iltemp, icol+2-ic,
     $                               c, s, a( jch-iskew*ic+ioffst, ic ),
     $                               ilda, extra, ctemp )
                        IF( iltemp ) THEN
                           CALL zlartg( a( jch-iskew*icol+ioffst,
     $                                  icol ), ctemp, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = realc*dummy
                           s = s*dummy
                           il = min( iendch, jch+jkl+jku ) + 2 - jch
                           extra = czero
                           CALL zlarot( .false., .true.,
     $                                  jch+jkl+jku.LE.iendch, il, c, s,
     $                                  a( jch-iskew*icol+ioffst,
     $                                  icol ), ilda, ctemp, extra )
                           ic = icol
                        END IF
  110                CONTINUE
  120             CONTINUE
  130          CONTINUE
*
               jku = uub
               DO 160 jkl = 1, llb
*
*                 Transform from bandwidth JKL-1, JKU to JKL, JKU
*
*                 First row actually rotated is MIN( N+JKL, M )
*                 First column actually rotated is N
*
                  iendch = min( n, m+jku ) - 1
                  DO 150 jr = min( n+jkl, m ) - 1, 1 - jku, -1
                     extra = czero
                     angle = twopi*dlarnd( 1, iseed )
                     c = cos( angle )*zlarnd( 5, iseed )
                     s = sin( angle )*zlarnd( 5, iseed )
                     icol = max( 1, jr-jkl+1 )
                     IF( jr.GT.0 ) THEN
                        il = min( n, jr+jku+1 ) + 1 - icol
                        CALL zlarot( .true., .false., jr+jku.LT.n, il,
     $                               c, s, a( jr-iskew*icol+ioffst,
     $                               icol ), ilda, dummy, extra )
                     END IF
*
*                    Chase "EXTRA" back down
*
                     ir = jr
                     DO 140 jch = jr + jku, iendch, jkl + jku
                        ilextr = ir.GT.0
                        IF( ilextr ) THEN
                           CALL zlartg( a( ir-iskew*jch+ioffst, jch ),
     $                                  extra, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = realc*dummy
                           s = s*dummy
                        END IF
                        ir = max( 1, ir )
                        irow = min( m-1, jch+jkl )
                        iltemp = jch + jkl.LT.m
                        ctemp = czero
                        CALL zlarot( .false., ilextr, iltemp, irow+2-ir,
     $                               c, s, a( ir-iskew*jch+ioffst,
     $                               jch ), ilda, extra, ctemp )
                        IF( iltemp ) THEN
                           CALL zlartg( a( irow-iskew*jch+ioffst, jch ),
     $                                  ctemp, realc, s, dummy )
                           dummy = zlarnd( 5, iseed )
                           c = realc*dummy
                           s = s*dummy
                           il = min( iendch, jch+jkl+jku ) + 2 - jch
                           extra = czero
                           CALL zlarot( .true., .true.,
     $                                  jch+jkl+jku.LE.iendch, il, c, s,
     $                                  a( irow-iskew*jch+ioffst, jch ),
     $                                  ilda, ctemp, extra )
                           ir = irow
                        END IF
  140                CONTINUE
  150             CONTINUE
  160          CONTINUE
*
            END IF
*
         ELSE
*
*           Symmetric -- A = U D U'
*           Hermitian -- A = U D U*
*
            ipackg = ipack
            ioffg = ioffst
*
            IF( topdwn ) THEN
*
*              Top-Down -- Generate Upper triangle only
*
               IF( ipack.GE.5 ) THEN
                  ipackg = 6
                  ioffg = uub + 1
               ELSE
                  ipackg = 1
               END IF
*
               DO 170 j = 1, mnmin
                  a( ( 1-iskew )*j+ioffg, j ) = dcmplx( d( j ) )
  170          CONTINUE
*
               DO 200 k = 1, uub
                  DO 190 jc = 1, n - 1
                     irow = max( 1, jc-k )
                     il = min( jc+1, k+2 )
                     extra = czero
                     ctemp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
                     angle = twopi*dlarnd( 1, iseed )
                     c = cos( angle )*zlarnd( 5, iseed )
                     s = sin( angle )*zlarnd( 5, iseed )
                     IF( zsym ) THEN
                        ct = c
                        st = s
                     ELSE
                        ctemp = dconjg( ctemp )
                        ct = dconjg( c )
                        st = dconjg( s )
                     END IF
                     CALL zlarot( .false., jc.GT.k, .true., il, c, s,
     $                            a( irow-iskew*jc+ioffg, jc ), ilda,
     $                            extra, ctemp )
                     CALL zlarot( .true., .true., .false.,
     $                            min( k, n-jc )+1, ct, st,
     $                            a( ( 1-iskew )*jc+ioffg, jc ), ilda,
     $                            ctemp, dummy )
*
*                    Chase EXTRA back up the matrix
*
                     icol = jc
                     DO 180 jch = jc - k, 1, -k
                        CALL zlartg( a( jch+1-iskew*( icol+1 )+ioffg,
     $                               icol+1 ), extra, realc, s, dummy )
                        dummy = zlarnd( 5, iseed )
                        c = dconjg( realc*dummy )
                        s = dconjg( -s*dummy )
                        ctemp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
                        IF( zsym ) THEN
                           ct = c
                           st = s
                        ELSE
                           ctemp = dconjg( ctemp )
                           ct = dconjg( c )
                           st = dconjg( s )
                        END IF
                        CALL zlarot( .true., .true., .true., k+2, c, s,
     $                               a( ( 1-iskew )*jch+ioffg, jch ),
     $                               ilda, ctemp, extra )
                        irow = max( 1, jch-k )
                        il = min( jch+1, k+2 )
                        extra = czero
                        CALL zlarot( .false., jch.GT.k, .true., il, ct,
     $                               st, a( irow-iskew*jch+ioffg, jch ),
     $                               ilda, extra, ctemp )
                        icol = jch
  180                CONTINUE
  190             CONTINUE
  200          CONTINUE
*
*              If we need lower triangle, copy from upper. Note that
*              the order of copying is chosen to work for 'q' -> 'b'
*
               IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
                  DO 230 jc = 1, n
                     irow = ioffst - iskew*jc
                     IF( zsym ) THEN
                        DO 210 jr = jc, min( n, jc+uub )
                           a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
  210                   CONTINUE
                     ELSE
                        DO 220 jr = jc, min( n, jc+uub )
                           a( jr+irow, jc ) = dconjg( a( jc-iskew*jr+
     $                                        ioffg, jr ) )
  220                   CONTINUE
                     END IF
  230             CONTINUE
                  IF( ipack.EQ.5 ) THEN
                     DO 250 jc = n - uub + 1, n
                        DO 240 jr = n + 2 - jc, uub + 1
                           a( jr, jc ) = czero
  240                   CONTINUE
  250                CONTINUE
                  END IF
                  IF( ipackg.EQ.6 ) THEN
                     ipackg = ipack
                  ELSE
                     ipackg = 0
                  END IF
               END IF
            ELSE
*
*              Bottom-Up -- Generate Lower triangle only
*
               IF( ipack.GE.5 ) THEN
                  ipackg = 5
                  IF( ipack.EQ.6 )
     $               ioffg = 1
               ELSE
                  ipackg = 2
               END IF
*
               DO 260 j = 1, mnmin
                  a( ( 1-iskew )*j+ioffg, j ) = dcmplx( d( j ) )
  260          CONTINUE
*
               DO 290 k = 1, uub
                  DO 280 jc = n - 1, 1, -1
                     il = min( n+1-jc, k+2 )
                     extra = czero
                     ctemp = a( 1+( 1-iskew )*jc+ioffg, jc )
                     angle = twopi*dlarnd( 1, iseed )
                     c = cos( angle )*zlarnd( 5, iseed )
                     s = sin( angle )*zlarnd( 5, iseed )
                     IF( zsym ) THEN
                        ct = c
                        st = s
                     ELSE
                        ctemp = dconjg( ctemp )
                        ct = dconjg( c )
                        st = dconjg( s )
                     END IF
                     CALL zlarot( .false., .true., n-jc.GT.k, il, c, s,
     $                            a( ( 1-iskew )*jc+ioffg, jc ), ilda,
     $                            ctemp, extra )
                     icol = max( 1, jc-k+1 )
                     CALL zlarot( .true., .false., .true., jc+2-icol,
     $                            ct, st, a( jc-iskew*icol+ioffg,
     $                            icol ), ilda, dummy, ctemp )
*
*                    Chase EXTRA back down the matrix
*
                     icol = jc
                     DO 270 jch = jc + k, n - 1, k
                        CALL zlartg( a( jch-iskew*icol+ioffg, icol ),
     $                               extra, realc, s, dummy )
                        dummy = zlarnd( 5, iseed )
                        c = realc*dummy
                        s = s*dummy
                        ctemp = a( 1+( 1-iskew )*jch+ioffg, jch )
                        IF( zsym ) THEN
                           ct = c
                           st = s
                        ELSE
                           ctemp = dconjg( ctemp )
                           ct = dconjg( c )
                           st = dconjg( s )
                        END IF
                        CALL zlarot( .true., .true., .true., k+2, c, s,
     $                               a( jch-iskew*icol+ioffg, icol ),
     $                               ilda, extra, ctemp )
                        il = min( n+1-jch, k+2 )
                        extra = czero
                        CALL zlarot( .false., .true., n-jch.GT.k, il,
     $                               ct, st, a( ( 1-iskew )*jch+ioffg,
     $                               jch ), ilda, ctemp, extra )
                        icol = jch
  270                CONTINUE
  280             CONTINUE
  290          CONTINUE
*
*              If we need upper triangle, copy from lower. Note that
*              the order of copying is chosen to work for 'b' -> 'q'
*
               IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
                  DO 320 jc = n, 1, -1
                     irow = ioffst - iskew*jc
                     IF( zsym ) THEN
                        DO 300 jr = jc, max( 1, jc-uub ), -1
                           a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
  300                   CONTINUE
                     ELSE
                        DO 310 jr = jc, max( 1, jc-uub ), -1
                           a( jr+irow, jc ) = dconjg( a( jc-iskew*jr+
     $                                        ioffg, jr ) )
  310                   CONTINUE
                     END IF
  320             CONTINUE
                  IF( ipack.EQ.6 ) THEN
                     DO 340 jc = 1, uub
                        DO 330 jr = 1, uub + 1 - jc
                           a( jr, jc ) = czero
  330                   CONTINUE
  340                CONTINUE
                  END IF
                  IF( ipackg.EQ.5 ) THEN
                     ipackg = ipack
                  ELSE
                     ipackg = 0
                  END IF
               END IF
            END IF
*
*           Ensure that the diagonal is real if Hermitian
*
            IF( .NOT.zsym ) THEN
               DO 350 jc = 1, n
                  irow = ioffst + ( 1-iskew )*jc
                  a( irow, jc ) = dcmplx( dble( a( irow, jc ) ) )
  350          CONTINUE
            END IF
*
         END IF
*
      ELSE
*
*        4)      Generate Banded Matrix by first
*                Rotating by random Unitary matrices,
*                then reducing the bandwidth using Householder
*                transformations.
*
*                Note: we should get here only if LDA .ge. N
*
         IF( isym.EQ.1 ) THEN
*
*           Non-symmetric -- A = U D V
*
            CALL zlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
     $                   iinfo )
         ELSE
*
*           Symmetric -- A = U D U' or
*           Hermitian -- A = U D U*
*
            IF( zsym ) THEN
               CALL zlagsy( m, llb, d, a, lda, iseed, work, iinfo )
            ELSE
               CALL zlaghe( m, llb, d, a, lda, iseed, work, iinfo )
            END IF
         END IF
*
         IF( iinfo.NE.0 ) THEN
            info = 3
            RETURN
         END IF
      END IF
*
*     5)      Pack the matrix
*
      IF( ipack.NE.ipackg ) THEN
         IF( ipack.EQ.1 ) THEN
*
*           'U' -- Upper triangular, not packed
*
            DO 370 j = 1, m
               DO 360 i = j + 1, m
                  a( i, j ) = czero
  360          CONTINUE
  370       CONTINUE
*
         ELSE IF( ipack.EQ.2 ) THEN
*
*           'L' -- Lower triangular, not packed
*
            DO 390 j = 2, m
               DO 380 i = 1, j - 1
                  a( i, j ) = czero
  380          CONTINUE
  390       CONTINUE
*
         ELSE IF( ipack.EQ.3 ) THEN
*
*           'C' -- Upper triangle packed Columnwise.
*
            icol = 1
            irow = 0
            DO 410 j = 1, m
               DO 400 i = 1, j
                  irow = irow + 1
                  IF( irow.GT.lda ) THEN
                     irow = 1
                     icol = icol + 1
                  END IF
                  a( irow, icol ) = a( i, j )
  400          CONTINUE
  410       CONTINUE
*
         ELSE IF( ipack.EQ.4 ) THEN
*
*           'R' -- Lower triangle packed Columnwise.
*
            icol = 1
            irow = 0
            DO 430 j = 1, m
               DO 420 i = j, m
                  irow = irow + 1
                  IF( irow.GT.lda ) THEN
                     irow = 1
                     icol = icol + 1
                  END IF
                  a( irow, icol ) = a( i, j )
  420          CONTINUE
  430       CONTINUE
*
         ELSE IF( ipack.GE.5 ) THEN
*
*           'B' -- The lower triangle is packed as a band matrix.
*           'Q' -- The upper triangle is packed as a band matrix.
*           'Z' -- The whole matrix is packed as a band matrix.
*
            IF( ipack.EQ.5 )
     $         uub = 0
            IF( ipack.EQ.6 )
     $         llb = 0
*
            DO 450 j = 1, uub
               DO 440 i = min( j+llb, m ), 1, -1
                  a( i-j+uub+1, j ) = a( i, j )
  440          CONTINUE
  450       CONTINUE
*
            DO 470 j = uub + 2, n
               DO 460 i = j - uub, min( j+llb, m )
                  a( i-j+uub+1, j ) = a( i, j )
  460          CONTINUE
  470       CONTINUE
         END IF
*
*        If packed, zero out extraneous elements.
*
*        Symmetric/Triangular Packed --
*        zero out everything after A(IROW,ICOL)
*
         IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
            DO 490 jc = icol, m
               DO 480 jr = irow + 1, lda
                  a( jr, jc ) = czero
  480          CONTINUE
               irow = 0
  490       CONTINUE
*
         ELSE IF( ipack.GE.5 ) THEN
*
*           Packed Band --
*              1st row is now in A( UUB+2-j, j), zero above it
*              m-th row is now in A( M+UUB-j,j), zero below it
*              last non-zero diagonal is now in A( UUB+LLB+1,j ),
*                 zero below it, too.
*
            ir1 = uub + llb + 2
            ir2 = uub + m + 2
            DO 520 jc = 1, n
               DO 500 jr = 1, uub + 1 - jc
                  a( jr, jc ) = czero
  500          CONTINUE
               DO 510 jr = max( 1, min( ir1, ir2-jc ) ), lda
                  a( jr, jc ) = czero
  510          CONTINUE
  520       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of ZLATMS
*
      END