15 files changed, 5706 insertions, 12 deletions

diff --git a/TESTING/MATGEN/Makefile b/TESTING/MATGEN/Makefile
index 350cb250..389b5540 100644
--- a/TESTING/MATGEN/Makefile
+++ b/TESTING/MATGEN/Makefile
@@ -35,23 +35,23 @@ include ../../make.inc
  
 SCATGEN = slatm1.o slaran.o slarnd.o
 
-SMATGEN = slatms.o slatme.o slatmr.o \
+SMATGEN = slatms.o slatme.o slatmr.o slatmt.o \
    slagge.o slagsy.o slakf2.o slarge.o slaror.o slarot.o slatm2.o \
-   slatm3.o slatm5.o slatm6.o
+   slatm3.o slatm5.o slatm6.o slatm7.o slahilb.o
 
-CMATGEN = clatms.o clatme.o clatmr.o \
+CMATGEN = clatms.o clatme.o clatmr.o clatmt.o \
    clagge.o claghe.o clagsy.o clakf2.o clarge.o claror.o clarot.o \
-   clatm1.o clarnd.o clatm2.o clatm3.o clatm5.o clatm6.o
+   clatm1.o clarnd.o clatm2.o clatm3.o clatm5.o clatm6.o clahilb.o
 
 DZATGEN = dlatm1.o dlaran.o dlarnd.o
 
-DMATGEN = dlatms.o dlatme.o dlatmr.o \
+DMATGEN = dlatms.o dlatme.o dlatmr.o dlatmt.o \
    dlagge.o dlagsy.o dlakf2.o dlarge.o dlaror.o dlarot.o dlatm2.o \
-   dlatm3.o dlatm5.o dlatm6.o
+   dlatm3.o dlatm5.o dlatm6.o dlatm7.o dlahilb.o
 
-ZMATGEN = zlatms.o zlatme.o zlatmr.o \
+ZMATGEN = zlatms.o zlatme.o zlatmr.o zlatmt.o \
    zlagge.o zlaghe.o zlagsy.o zlakf2.o zlarge.o zlaror.o zlarot.o \
-   zlatm1.o zlarnd.o zlatm2.o zlatm3.o zlatm5.o zlatm6.o
+   zlatm1.o zlarnd.o zlatm2.o zlatm3.o zlatm5.o zlatm6.o zlahilb.o
 
 all:	../../$(TMGLIB)
 
diff --git a/TESTING/MATGEN/clahilb.f b/TESTING/MATGEN/clahilb.f
new file mode 100644
index 00000000..fcc054eb
--- /dev/null
+++ b/TESTING/MATGEN/clahilb.f
@@ -0,0 +1,210 @@
+      SUBROUTINE CLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, 
+     $     INFO, PATH)
+!
+!  -- LAPACK auxiliary test routine (version 3.0) --
+!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+!     Courant Institute, Argonne National Lab, and Rice University
+!     28 August, 2006
+!
+!     David Vu <dtv@cs.berkeley.edu>
+!     Yozo Hida <yozo@cs.berkeley.edu>
+!     Jason Riedy <ejr@cs.berkeley.edu>
+!     D. Halligan <dhalligan@berkeley.edu>
+!
+      IMPLICIT NONE
+!     .. Scalar Arguments ..
+      INTEGER N, NRHS, LDA, LDX, LDB, INFO
+!     .. Array Arguments ..
+      REAL WORK(N)
+      COMPLEX A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
+      CHARACTER*3        PATH
+!     ..
+!
+!  Purpose
+!  =======
+!
+!  CLAHILB generates an N by N scaled Hilbert matrix in A along with
+!  NRHS right-hand sides in B and solutions in X such that A*X=B.
+!
+!  The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
+!  entries are integers.  The right-hand sides are the first NRHS
+!  columns of M * the identity matrix, and the solutions are the
+!  first NRHS columns of the inverse Hilbert matrix.
+!
+!  The condition number of the Hilbert matrix grows exponentially with
+!  its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
+!  Hilbert matrices beyond a relatively small dimension cannot be
+!  generated exactly without extra precision.  Precision is exhausted
+!  when the largest entry in the inverse Hilbert matrix is greater than
+!  2 to the power of the number of bits in the fraction of the data type
+!  used plus one, which is 24 for single precision.
+!
+!  In single, the generated solution is exact for N <= 6 and has
+!  small componentwise error for 7 <= N <= 11.
+!
+!  Arguments
+!  =========
+!
+!  N       (input) INTEGER
+!          The dimension of the matrix A.
+!
+!  NRHS    (input) NRHS
+!          The requested number of right-hand sides.
+!
+!  A       (output) COMPLEX array, dimension (LDA, N)
+!          The generated scaled Hilbert matrix.
+!
+!  LDA     (input) INTEGER
+!          The leading dimension of the array A.  LDA >= N.
+!
+!  X       (output) COMPLEX array, dimension (LDX, NRHS)
+!          The generated exact solutions.  Currently, the first NRHS
+!          columns of the inverse Hilbert matrix.
+!
+!  LDX     (input) INTEGER
+!          The leading dimension of the array X.  LDX >= N.
+!
+!  B       (output) REAL array, dimension (LDB, NRHS)
+!          The generated right-hand sides.  Currently, the first NRHS
+!          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
+!
+!  LDB     (input) INTEGER
+!          The leading dimension of the array B.  LDB >= N.
+!
+!  WORK    (workspace) REAL array, dimension (N)
+!
+!
+!  INFO    (output) INTEGER
+!          = 0: successful exit
+!          = 1: N is too large; the data is still generated but may not
+!               be not exact.
+!          < 0: if INFO = -i, the i-th argument had an illegal value
+!
+!  =====================================================================
+
+!     .. Local Scalars ..
+      INTEGER TM, TI, R
+      INTEGER M
+      INTEGER I, J
+      COMPLEX TMP
+      CHARACTER*2 C2
+
+!     .. Parameters ..
+!     NMAX_EXACT   the largest dimension where the generated data is
+!                  exact.
+!     NMAX_APPROX  the largest dimension where the generated data has
+!                  a small componentwise relative error.
+!     ??? complex uses how many bits ???
+      INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
+      PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8)
+
+!     d's are generated from random permuation of those eight elements.
+      COMPLEX D1(8), D2(8), INVD1(8), INVD2(8) 
+      DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
+      DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
+      
+      DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
+     $     (-.5,-.5),(.5,-.5),(.5,.5)/
+      DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
+     $     (-.5,.5),(.5,.5),(.5,-.5)/
+      
+!     ..
+!     .. External Functions
+      EXTERNAL CLASET, LSAMEN
+      INTRINSIC REAL
+      LOGICAL LSAMEN
+!     ..
+!     .. Executable Statements ..
+      C2 = PATH( 2: 3 )
+!
+!     Test the input arguments
+!
+      INFO = 0
+      IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
+         INFO = -1
+      ELSE IF (NRHS .LT. 0) THEN
+         INFO = -2
+      ELSE IF (LDA .LT. N) THEN
+         INFO = -4
+      ELSE IF (LDX .LT. N) THEN
+         INFO = -6
+      ELSE IF (LDB .LT. N) THEN
+         INFO = -8
+      END IF
+      IF (INFO .LT. 0) THEN
+         CALL XERBLA('CLAHILB', -INFO)
+         RETURN
+      END IF
+      IF (N .GT. NMAX_EXACT) THEN
+         INFO = 1
+      END IF
+
+!     Compute M = the LCM of the integers [1, 2*N-1].  The largest
+!     reasonable N is small enough that integers suffice (up to N = 11).
+      M = 1
+      DO I = 2, (2*N-1)
+         TM = M
+         TI = I
+         R = MOD(TM, TI)
+         DO WHILE (R .NE. 0)
+            TM = TI
+            TI = R
+            R = MOD(TM, TI)
+         END DO
+         M = (M / TI) * I
+      END DO
+
+!     Generate the scaled Hilbert matrix in A
+!     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i*
+      IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
+         DO J = 1, N
+            DO I = 1, N
+               A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
+     $              * D1(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      ELSE
+         DO J = 1, N
+            DO I = 1, N
+               A(I, J) = D1(MOD(J,SIZE_D)+1) * (REAL(M) / (I + J - 1))
+     $              * D2(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      END IF
+
+!     Generate matrix B as simply the first NRHS columns of M * the
+!     identity.
+      TMP = REAL(M)
+      CALL CLASET('Full', N, NRHS, (0.0,0.0), TMP, B, LDB)
+
+!     Generate the true solutions in X.  Because B = the first NRHS
+!     columns of M*I, the true solutions are just the first NRHS columns
+!     of the inverse Hilbert matrix.
+      WORK(1) = N
+      DO J = 2, N
+         WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
+     $        * (N +J -1)
+      END DO
+      
+!     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i*
+      IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
+         DO J = 1, NRHS
+            DO I = 1, N
+               X(I, J) = 
+     $              INVD1(MOD(J,SIZE_D)+1) *
+     $              ((WORK(I)*WORK(J)) / (I + J - 1)) 
+     $              * INVD1(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      ELSE
+         DO J = 1, NRHS
+            DO I = 1, N
+               X(I, J) =
+     $              INVD2(MOD(J,SIZE_D)+1) *
+     $              ((WORK(I)*WORK(J)) / (I + J - 1))
+     $              * INVD1(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      END IF
+      END
+      
diff --git a/TESTING/MATGEN/clarot.f b/TESTING/MATGEN/clarot.f
index bdf0a009..f73d10b4 100644
--- a/TESTING/MATGEN/clarot.f
+++ b/TESTING/MATGEN/clarot.f
@@ -19,7 +19,6 @@
 *
 *     CLAROT applies a (Givens) rotation to two adjacent rows or
 *     columns, where one element of the first and/or last column/row
-*     November 2006
 *     for use on matrices stored in some format other than GE, so
 *     that elements of the matrix may be used or modified for which
 *     no array element is provided.
diff --git a/TESTING/MATGEN/clatmt.f b/TESTING/MATGEN/clatmt.f
new file mode 100644
index 00000000..a1b8e372
--- /dev/null
+++ b/TESTING/MATGEN/clatmt.f
@@ -0,0 +1,1172 @@
+      SUBROUTINE CLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
+     $                   RANK, KL, KU, PACK, A, LDA, WORK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Craig Lucas, University of Manchester / NAG Ltd.
+*     October, 2008
+*
+*     .. Scalar Arguments ..
+      REAL               COND, DMAX
+      INTEGER            INFO, KL, KU, LDA, M, MODE, N, RANK
+      CHARACTER          DIST, PACK, SYM
+*     ..
+*     .. Array Arguments ..
+      COMPLEX            A( LDA, * ), WORK( * )
+      REAL               D( * )
+      INTEGER            ISEED( 4 )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*     CLATMT generates random matrices with specified singular values
+*     (or hermitian with specified eigenvalues)
+*     for testing LAPACK programs.
+*
+*     CLATMT 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 CLATMT
+*           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      - REAL array, dimension ( MIN( M, N ) )
+*           This array is used to specify the singular values or
+*           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
+*           assumed to contain the singular/eigenvalues, otherwise
+*           they will be computed according to MODE, COND, and DMAX,
+*           and placed in D.
+*           Modified if MODE is nonzero.
+*
+*  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:RANK)=1.0/COND
+*           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
+*           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
+*           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
+*           MODE = 5 sets D to random numbers in the range
+*                    ( 1/COND , 1 ) such that their logarithms
+*                    are uniformly distributed.
+*           MODE = 6 set D to random numbers from same distribution
+*                    as the rest of the matrix.
+*           MODE < 0 has the same meaning as ABS(MODE), except that
+*              the order of the elements of D is reversed.
+*           Thus if MODE is positive, D has entries ranging from
+*              1 to 1/COND, if negative, from 1/COND to 1,
+*           If SYM='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   - REAL
+*           On entry, this is used as described under MODE above.
+*           If used, it must be >= 1. Not modified.
+*
+*  DMAX   - REAL
+*           If MODE is neither -6, 0 nor 6, the contents of D, as
+*           computed according to MODE and COND, will be scaled by
+*           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
+*           singular value (which is to say the norm) will be abs(DMAX).
+*           Note that DMAX need not be positive: if DMAX is negative
+*           (or zero), D will be scaled by a negative number (or zero).
+*           Not modified.
+*
+*  RANK   - INTEGER
+*           The rank of matrix to be generated for modes 1,2,3 only.
+*           D( RANK+1:N ) = 0.
+*           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 CLATMT differ only in the PACK parameter,
+*           they will generate mathematically equivalent matrices.
+*           Not modified.
+*
+*  A      - COMPLEX 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 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 SLATM7
+*            2  => Cannot scale to DMAX (max. sing. value is 0)
+*            3  => Error return from CLAGGE, CLAGHE or CLAGSY
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E+0 )
+      REAL               ONE
+      PARAMETER          ( ONE = 1.0E+0 )
+      COMPLEX            CZERO
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
+      REAL               TWOPI
+      PARAMETER          ( TWOPI = 6.2831853071795864769252867663E+0 )
+*     ..
+*     .. Local Scalars ..
+      COMPLEX            C, CT, CTEMP, DUMMY, EXTRA, S, ST
+      REAL               ALPHA, ANGLE, REALC, TEMP
+      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
+      LOGICAL            CSYM, GIVENS, ILEXTR, ILTEMP, TOPDWN
+*     ..
+*     .. External Functions ..
+      COMPLEX            CLARND
+      REAL               SLARND
+      LOGICAL            LSAME
+      EXTERNAL           CLARND, SLARND, LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLAGGE, CLAGHE, CLAGSY, CLAROT, CLARTG, CLASET,
+     $                   SLATM7, SSCAL, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, CMPLX, CONJG, COS, MAX, MIN, MOD, REAL,
+     $                   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
+         CSYM = .FALSE.
+      ELSE IF( LSAME( SYM, 'P' ) ) THEN
+         ISYM = 2
+         IRSIGN = 0
+         CSYM = .FALSE.
+      ELSE IF( LSAME( SYM, 'S' ) ) THEN
+         ISYM = 2
+         IRSIGN = 0
+         CSYM = .TRUE.
+      ELSE IF( LSAME( SYM, 'H' ) ) THEN
+         ISYM = 2
+         IRSIGN = 1
+         CSYM = .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 )
+*
+      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( REAL( LLB+UUB ).LT.0.3*REAL( 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( 'CLATMT', -INFO )
+         RETURN
+      END IF
+*
+*     Initialize random number generator
+*
+      DO 100 I = 1, 4
+         ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
+  100 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 SLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, RANK,
+     $             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( RANK ) ) ) 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 110 I = 2, RANK
+            TEMP = MAX( TEMP, ABS( D( I ) ) )
+  110    CONTINUE
+*
+         IF( TEMP.GT.ZERO ) THEN
+            ALPHA = DMAX / TEMP
+         ELSE
+            INFO = 2
+            RETURN
+         END IF
+*
+         CALL SSCAL( RANK, ALPHA, D, 1 )
+*
+      END IF
+*
+      CALL CLASET( '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 120 J = 1, MNMIN
+            A( ( 1-ISKEW )*J+IOFFST, J ) = CMPLX( D( J ) )
+  120    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 130 J = 1, MNMIN
+               A( ( 1-ISKEW )*J+IOFFST, J ) = CMPLX( D( J ) )
+  130       CONTINUE
+*
+            IF( TOPDWN ) THEN
+               JKL = 0
+               DO 160 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 150 JR = 1, MIN( M+JKU, N ) + JKL - 1
+                     EXTRA = CZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )*CLARND( 5, ISEED )
+                     S = SIN( ANGLE )*CLARND( 5, ISEED )
+                     ICOL = MAX( 1, JR-JKL )
+                     IF( JR.LT.M ) THEN
+                        IL = MIN( N, JR+JKU ) + 1 - ICOL
+                        CALL CLAROT( .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 140 JCH = JR - JKL, 1, -JKL - JKU
+                        IF( IR.LT.M ) THEN
+                           CALL CLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), EXTRA, REALC, S, DUMMY )
+                           DUMMY = CLARND( 5, ISEED )
+                           C = CONJG( REALC*DUMMY )
+                           S = CONJG( -S*DUMMY )
+                        END IF
+                        IROW = MAX( 1, JCH-JKU )
+                        IL = IR + 2 - IROW
+                        CTEMP = CZERO
+                        ILTEMP = JCH.GT.JKU
+                        CALL CLAROT( .FALSE., ILTEMP, .TRUE., IL, C, S,
+     $                               A( IROW-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, CTEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL CLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), CTEMP, REALC, S, DUMMY )
+                           DUMMY = CLARND( 5, ISEED )
+                           C = CONJG( REALC*DUMMY )
+                           S = CONJG( -S*DUMMY )
+*
+                           ICOL = MAX( 1, JCH-JKU-JKL )
+                           IL = IC + 2 - ICOL
+                           EXTRA = CZERO
+                           CALL CLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
+     $                                  IL, C, S, A( IROW-ISKEW*ICOL+
+     $                                  IOFFST, ICOL ), ILDA, EXTRA,
+     $                                  CTEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  140                CONTINUE
+  150             CONTINUE
+  160          CONTINUE
+*
+               JKU = UUB
+               DO 190 JKL = 1, LLB
+*
+*                 Transform from bandwidth JKL-1, JKU to JKL, JKU
+*
+                  DO 180 JC = 1, MIN( N+JKL, M ) + JKU - 1
+                     EXTRA = CZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )*CLARND( 5, ISEED )
+                     S = SIN( ANGLE )*CLARND( 5, ISEED )
+                     IROW = MAX( 1, JC-JKU )
+                     IF( JC.LT.N ) THEN
+                        IL = MIN( M, JC+JKL ) + 1 - IROW
+                        CALL CLAROT( .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 170 JCH = JC - JKU, 1, -JKL - JKU
+                        IF( IC.LT.N ) THEN
+                           CALL CLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), EXTRA, REALC, S, DUMMY )
+                           DUMMY = CLARND( 5, ISEED )
+                           C = CONJG( REALC*DUMMY )
+                           S = CONJG( -S*DUMMY )
+                        END IF
+                        ICOL = MAX( 1, JCH-JKL )
+                        IL = IC + 2 - ICOL
+                        CTEMP = CZERO
+                        ILTEMP = JCH.GT.JKL
+                        CALL CLAROT( .TRUE., ILTEMP, .TRUE., IL, C, S,
+     $                               A( IR-ISKEW*ICOL+IOFFST, ICOL ),
+     $                               ILDA, CTEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL CLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
+     $                                  ICOL+1 ), CTEMP, REALC, S,
+     $                                  DUMMY )
+                           DUMMY = CLARND( 5, ISEED )
+                           C = CONJG( REALC*DUMMY )
+                           S = CONJG( -S*DUMMY )
+                           IROW = MAX( 1, JCH-JKL-JKU )
+                           IL = IR + 2 - IROW
+                           EXTRA = CZERO
+                           CALL CLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
+     $                                  IL, C, S, A( IROW-ISKEW*ICOL+
+     $                                  IOFFST, ICOL ), ILDA, EXTRA,
+     $                                  CTEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  170                CONTINUE
+  180             CONTINUE
+  190          CONTINUE
+*
+            ELSE
+*
+*              Bottom-Up -- Start at the bottom right.
+*
+               JKL = 0
+               DO 220 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 210 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
+                     EXTRA = CZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )*CLARND( 5, ISEED )
+                     S = SIN( ANGLE )*CLARND( 5, ISEED )
+                     IROW = MAX( 1, JC-JKU+1 )
+                     IF( JC.GT.0 ) THEN
+                        IL = MIN( M, JC+JKL+1 ) + 1 - IROW
+                        CALL CLAROT( .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 200 JCH = JC + JKL, IENDCH, JKL + JKU
+                        ILEXTR = IC.GT.0
+                        IF( ILEXTR ) THEN
+                           CALL CLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                                  EXTRA, REALC, S, DUMMY )
+                           DUMMY = CLARND( 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 CLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
+     $                               C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, EXTRA, CTEMP )
+                        IF( ILTEMP ) THEN
+                           CALL CLARTG( A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), CTEMP, REALC, S, DUMMY )
+                           DUMMY = CLARND( 5, ISEED )
+                           C = REALC*DUMMY
+                           S = S*DUMMY
+                           IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+                           EXTRA = CZERO
+                           CALL CLAROT( .FALSE., .TRUE.,
+     $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+     $                                  A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), ILDA, CTEMP, EXTRA )
+                           IC = ICOL
+                        END IF
+  200                CONTINUE
+  210             CONTINUE
+  220          CONTINUE
+*
+               JKU = UUB
+               DO 250 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 240 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
+                     EXTRA = CZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )*CLARND( 5, ISEED )
+                     S = SIN( ANGLE )*CLARND( 5, ISEED )
+                     ICOL = MAX( 1, JR-JKL+1 )
+                     IF( JR.GT.0 ) THEN
+                        IL = MIN( N, JR+JKU+1 ) + 1 - ICOL
+                        CALL CLAROT( .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 230 JCH = JR + JKU, IENDCH, JKL + JKU
+                        ILEXTR = IR.GT.0
+                        IF( ILEXTR ) THEN
+                           CALL CLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ),
+     $                                  EXTRA, REALC, S, DUMMY )
+                           DUMMY = CLARND( 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 CLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
+     $                               C, S, A( IR-ISKEW*JCH+IOFFST,
+     $                               JCH ), ILDA, EXTRA, CTEMP )
+                        IF( ILTEMP ) THEN
+                           CALL CLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  CTEMP, REALC, S, DUMMY )
+                           DUMMY = CLARND( 5, ISEED )
+                           C = REALC*DUMMY
+                           S = S*DUMMY
+                           IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+                           EXTRA = CZERO
+                           CALL CLAROT( .TRUE., .TRUE.,
+     $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+     $                                  A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  ILDA, CTEMP, EXTRA )
+                           IR = IROW
+                        END IF
+  230                CONTINUE
+  240             CONTINUE
+  250          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 260 J = 1, MNMIN
+                  A( ( 1-ISKEW )*J+IOFFG, J ) = CMPLX( D( J ) )
+  260          CONTINUE
+*
+               DO 290 K = 1, UUB
+                  DO 280 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*SLARND( 1, ISEED )
+                     C = COS( ANGLE )*CLARND( 5, ISEED )
+                     S = SIN( ANGLE )*CLARND( 5, ISEED )
+                     IF( CSYM ) THEN
+                        CT = C
+                        ST = S
+                     ELSE
+                        CTEMP = CONJG( CTEMP )
+                        CT = CONJG( C )
+                        ST = CONJG( S )
+                     END IF
+                     CALL CLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
+     $                            A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
+     $                            EXTRA, CTEMP )
+                     CALL CLAROT( .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 270 JCH = JC - K, 1, -K
+                        CALL CLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
+     $                               ICOL+1 ), EXTRA, REALC, S, DUMMY )
+                        DUMMY = CLARND( 5, ISEED )
+                        C = CONJG( REALC*DUMMY )
+                        S = CONJG( -S*DUMMY )
+                        CTEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
+                        IF( CSYM ) THEN
+                           CT = C
+                           ST = S
+                        ELSE
+                           CTEMP = CONJG( CTEMP )
+                           CT = CONJG( C )
+                           ST = CONJG( S )
+                        END IF
+                        CALL CLAROT( .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 CLAROT( .FALSE., JCH.GT.K, .TRUE., IL, CT,
+     $                               ST, A( IROW-ISKEW*JCH+IOFFG, JCH ),
+     $                               ILDA, EXTRA, CTEMP )
+                        ICOL = JCH
+  270                CONTINUE
+  280             CONTINUE
+  290          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 320 JC = 1, N
+                     IROW = IOFFST - ISKEW*JC
+                     IF( CSYM ) THEN
+                        DO 300 JR = JC, MIN( N, JC+UUB )
+                           A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  300                   CONTINUE
+                     ELSE
+                        DO 310 JR = JC, MIN( N, JC+UUB )
+                           A( JR+IROW, JC ) = CONJG( A( JC-ISKEW*JR+
+     $                                        IOFFG, JR ) )
+  310                   CONTINUE
+                     END IF
+  320             CONTINUE
+                  IF( IPACK.EQ.5 ) THEN
+                     DO 340 JC = N - UUB + 1, N
+                        DO 330 JR = N + 2 - JC, UUB + 1
+                           A( JR, JC ) = CZERO
+  330                   CONTINUE
+  340                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 350 J = 1, MNMIN
+                  A( ( 1-ISKEW )*J+IOFFG, J ) = CMPLX( D( J ) )
+  350          CONTINUE
+*
+               DO 380 K = 1, UUB
+                  DO 370 JC = N - 1, 1, -1
+                     IL = MIN( N+1-JC, K+2 )
+                     EXTRA = CZERO
+                     CTEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )*CLARND( 5, ISEED )
+                     S = SIN( ANGLE )*CLARND( 5, ISEED )
+                     IF( CSYM ) THEN
+                        CT = C
+                        ST = S
+                     ELSE
+                        CTEMP = CONJG( CTEMP )
+                        CT = CONJG( C )
+                        ST = CONJG( S )
+                     END IF
+                     CALL CLAROT( .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 CLAROT( .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 360 JCH = JC + K, N - 1, K
+                        CALL CLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+     $                               EXTRA, REALC, S, DUMMY )
+                        DUMMY = CLARND( 5, ISEED )
+                        C = REALC*DUMMY
+                        S = S*DUMMY
+                        CTEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
+                        IF( CSYM ) THEN
+                           CT = C
+                           ST = S
+                        ELSE
+                           CTEMP = CONJG( CTEMP )
+                           CT = CONJG( C )
+                           ST = CONJG( S )
+                        END IF
+                        CALL CLAROT( .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 CLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL,
+     $                               CT, ST, A( ( 1-ISKEW )*JCH+IOFFG,
+     $                               JCH ), ILDA, CTEMP, EXTRA )
+                        ICOL = JCH
+  360                CONTINUE
+  370             CONTINUE
+  380          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 410 JC = N, 1, -1
+                     IROW = IOFFST - ISKEW*JC
+                     IF( CSYM ) THEN
+                        DO 390 JR = JC, MAX( 1, JC-UUB ), -1
+                           A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  390                   CONTINUE
+                     ELSE
+                        DO 400 JR = JC, MAX( 1, JC-UUB ), -1
+                           A( JR+IROW, JC ) = CONJG( A( JC-ISKEW*JR+
+     $                                        IOFFG, JR ) )
+  400                   CONTINUE
+                     END IF
+  410             CONTINUE
+                  IF( IPACK.EQ.6 ) THEN
+                     DO 430 JC = 1, UUB
+                        DO 420 JR = 1, UUB + 1 - JC
+                           A( JR, JC ) = CZERO
+  420                   CONTINUE
+  430                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.CSYM ) THEN
+               DO 440 JC = 1, N
+                  IROW = IOFFST + ( 1-ISKEW )*JC
+                  A( IROW, JC ) = CMPLX( REAL( A( IROW, JC ) ) )
+  440          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 CLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
+     $                   IINFO )
+         ELSE
+*
+*           Symmetric -- A = U D U' or
+*           Hermitian -- A = U D U*
+*
+            IF( CSYM ) THEN
+               CALL CLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
+            ELSE
+               CALL CLAGHE( 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 460 J = 1, M
+               DO 450 I = J + 1, M
+                  A( I, J ) = CZERO
+  450          CONTINUE
+  460       CONTINUE
+*
+         ELSE IF( IPACK.EQ.2 ) THEN
+*
+*           'L' -- Lower triangular, not packed
+*
+            DO 480 J = 2, M
+               DO 470 I = 1, J - 1
+                  A( I, J ) = CZERO
+  470          CONTINUE
+  480       CONTINUE
+*
+         ELSE IF( IPACK.EQ.3 ) THEN
+*
+*           'C' -- Upper triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 500 J = 1, M
+               DO 490 I = 1, J
+                  IROW = IROW + 1
+                  IF( IROW.GT.LDA ) THEN
+                     IROW = 1
+                     ICOL = ICOL + 1
+                  END IF
+                  A( IROW, ICOL ) = A( I, J )
+  490          CONTINUE
+  500       CONTINUE
+*
+         ELSE IF( IPACK.EQ.4 ) THEN
+*
+*           'R' -- Lower triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 520 J = 1, M
+               DO 510 I = J, M
+                  IROW = IROW + 1
+                  IF( IROW.GT.LDA ) THEN
+                     IROW = 1
+                     ICOL = ICOL + 1
+                  END IF
+                  A( IROW, ICOL ) = A( I, J )
+  510          CONTINUE
+  520       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 540 J = 1, UUB
+               DO 530 I = MIN( J+LLB, M ), 1, -1
+                  A( I-J+UUB+1, J ) = A( I, J )
+  530          CONTINUE
+  540       CONTINUE
+*
+            DO 560 J = UUB + 2, N
+               DO 550 I = J - UUB, MIN( J+LLB, M )
+                  A( I-J+UUB+1, J ) = A( I, J )
+  550          CONTINUE
+  560       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 580 JC = ICOL, M
+               DO 570 JR = IROW + 1, LDA
+                  A( JR, JC ) = CZERO
+  570          CONTINUE
+               IROW = 0
+  580       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 610 JC = 1, N
+               DO 590 JR = 1, UUB + 1 - JC
+                  A( JR, JC ) = CZERO
+  590          CONTINUE
+               DO 600 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
+                  A( JR, JC ) = CZERO
+  600          CONTINUE
+  610       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of CLATMT
+*
+      END
diff --git a/TESTING/MATGEN/dlahilb.f b/TESTING/MATGEN/dlahilb.f
new file mode 100644
index 00000000..ebc4d55b
--- /dev/null
+++ b/TESTING/MATGEN/dlahilb.f
@@ -0,0 +1,168 @@
+      SUBROUTINE DLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
+!
+!  -- LAPACK auxiliary test routine (version 3.0) --
+!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+!     Courant Institute, Argonne National Lab, and Rice University
+!     28 August, 2006
+!
+!     David Vu <dtv@cs.berkeley.edu>      
+!     Yozo Hida <yozo@cs.berkeley.edu>      
+!     Jason Riedy <ejr@cs.berkeley.edu>
+!     D. Halligan <dhalligan@berkeley.edu>
+!
+      IMPLICIT NONE
+!     .. Scalar Arguments ..
+      INTEGER N, NRHS, LDA, LDX, LDB, INFO
+!     .. Array Arguments ..
+      DOUBLE PRECISION A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
+!     ..
+!
+!  Purpose
+!  =======
+!
+!  DLAHILB generates an N by N scaled Hilbert matrix in A along with
+!  NRHS right-hand sides in B and solutions in X such that A*X=B.
+!
+!  The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
+!  entries are integers.  The right-hand sides are the first NRHS 
+!  columns of M * the identity matrix, and the solutions are the 
+!  first NRHS columns of the inverse Hilbert matrix.
+!
+!  The condition number of the Hilbert matrix grows exponentially with
+!  its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
+!  Hilbert matrices beyond a relatively small dimension cannot be
+!  generated exactly without extra precision.  Precision is exhausted
+!  when the largest entry in the inverse Hilbert matrix is greater than
+!  2 to the power of the number of bits in the fraction of the data type
+!  used plus one, which is 24 for single precision.  
+!
+!  In single, the generated solution is exact for N <= 6 and has
+!  small componentwise error for 7 <= N <= 11.
+!
+!  Arguments
+!  =========
+!
+!  N       (input) INTEGER
+!          The dimension of the matrix A.
+!      
+!  NRHS    (input) NRHS
+!          The requested number of right-hand sides.
+!
+!  A       (output) DOUBLE PRECISION array, dimension (LDA, N)
+!          The generated scaled Hilbert matrix.
+!
+!  LDA     (input) INTEGER
+!          The leading dimension of the array A.  LDA >= N.
+!
+!  X       (output) DOUBLE PRECISION array, dimension (LDX, NRHS)
+!          The generated exact solutions.  Currently, the first NRHS
+!          columns of the inverse Hilbert matrix.
+!
+!  LDX     (input) INTEGER
+!          The leading dimension of the array X.  LDX >= N.
+!
+!  B       (output) DOUBLE PRECISION array, dimension (LDB, NRHS)
+!          The generated right-hand sides.  Currently, the first NRHS
+!          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
+!
+!  LDB     (input) INTEGER
+!          The leading dimension of the array B.  LDB >= N.
+!
+!  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
+!
+!
+!  INFO    (output) INTEGER
+!          = 0: successful exit
+!          = 1: N is too large; the data is still generated but may not
+!               be not exact.
+!          < 0: if INFO = -i, the i-th argument had an illegal value
+!
+!  =====================================================================
+
+!     .. Local Scalars ..
+      INTEGER TM, TI, R
+      INTEGER M
+      INTEGER I, J
+      COMPLEX*16 TMP
+
+!     .. Parameters ..
+!     NMAX_EXACT   the largest dimension where the generated data is
+!                  exact.
+!     NMAX_APPROX  the largest dimension where the generated data has
+!                  a small componentwise relative error.
+      INTEGER NMAX_EXACT, NMAX_APPROX
+      PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
+
+!     ..
+!     .. External Functions
+      EXTERNAL DLASET
+      INTRINSIC DBLE
+!     ..
+!     .. Executable Statements ..
+!
+!     Test the input arguments
+!
+      INFO = 0
+      IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
+         INFO = -1
+      ELSE IF (NRHS .LT. 0) THEN
+         INFO = -2
+      ELSE IF (LDA .LT. N) THEN
+         INFO = -4
+      ELSE IF (LDX .LT. N) THEN
+         INFO = -6
+      ELSE IF (LDB .LT. N) THEN
+         INFO = -8
+      END IF
+      IF (INFO .LT. 0) THEN
+         CALL XERBLA('DLAHILB', -INFO)
+         RETURN
+      END IF
+      IF (N .GT. NMAX_EXACT) THEN
+         INFO = 1
+      END IF
+
+!     Compute M = the LCM of the integers [1, 2*N-1].  The largest
+!     reasonable N is small enough that integers suffice (up to N = 11).
+      M = 1
+      DO I = 2, (2*N-1)
+         TM = M
+         TI = I
+         R = MOD(TM, TI)
+         DO WHILE (R .NE. 0)
+            TM = TI
+            TI = R
+            R = MOD(TM, TI)
+         END DO
+         M = (M / TI) * I
+      END DO
+
+!     Generate the scaled Hilbert matrix in A
+      DO J = 1, N
+         DO I = 1, N
+            A(I, J) = DBLE(M) / (I + J - 1)
+         END DO
+      END DO
+
+!     Generate matrix B as simply the first NRHS columns of M * the
+!     identity.
+      TMP = DBLE(M)
+      CALL DLASET('Full', N, NRHS, 0.0D+0, TMP, B, LDB)
+
+!     Generate the true solutions in X.  Because B = the first NRHS
+!     columns of M*I, the true solutions are just the first NRHS columns
+!     of the inverse Hilbert matrix.
+      WORK(1) = N
+      DO J = 2, N
+         WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
+     $        * (N +J -1)
+      END DO
+      
+      DO J = 1, NRHS
+         DO I = 1, N
+            X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
+         END DO
+      END DO
+
+      END
+
diff --git a/TESTING/MATGEN/dlarot.f b/TESTING/MATGEN/dlarot.f
index 89842297..319ab1a8 100644
--- a/TESTING/MATGEN/dlarot.f
+++ b/TESTING/MATGEN/dlarot.f
@@ -19,7 +19,6 @@
 *
 *     DLAROT applies a (Givens) rotation to two adjacent rows or
 *     columns, where one element of the first and/or last column/row
-*     November 2006
 *     for use on matrices stored in some format other than GE, so
 *     that elements of the matrix may be used or modified for which
 *     no array element is provided.
diff --git a/TESTING/MATGEN/dlatm7.f b/TESTING/MATGEN/dlatm7.f
new file mode 100644
index 00000000..e21ad3c9
--- /dev/null
+++ b/TESTING/MATGEN/dlatm7.f
@@ -0,0 +1,255 @@
+      SUBROUTINE DLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, N,
+     $                   RANK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Craig Lucas, University of Manchester / NAG Ltd.
+*     October, 2008
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   COND
+      INTEGER            IDIST, INFO, IRSIGN, MODE, N, RANK
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   D( * )
+      INTEGER            ISEED( 4 )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*     DLATM7 computes the entries of D as specified by MODE
+*     COND and IRSIGN. IDIST and ISEED determine the generation
+*     of random numbers. DLATM7 is called by DLATMT to generate
+*     random test matrices.
+*
+*  Arguments
+*  =========
+*
+*  MODE   - INTEGER
+*           On entry describes how D is to be computed:
+*           MODE = 0 means do not change D.
+*
+*           MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
+*           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
+*           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1)) I=1:RANK
+*
+*           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,
+*           Not modified.
+*
+*  COND   - DOUBLE PRECISION
+*           On entry, used as described under MODE above.
+*           If used, it must be >= 1. Not modified.
+*
+*  IRSIGN - INTEGER
+*           On entry, if MODE neither -6, 0 nor 6, determines sign of
+*           entries of D
+*           0 => leave entries of D unchanged
+*           1 => multiply each entry of D by 1 or -1 with probability .5
+*
+*  IDIST  - CHARACTER*1
+*           On entry, IDIST specifies the type of distribution to be
+*           used to generate a random matrix .
+*           1 => UNIFORM( 0, 1 )
+*           2 => UNIFORM( -1, 1 )
+*           3 => NORMAL( 0, 1 )
+*           Not modified.
+*
+*  ISEED  - INTEGER array, dimension ( 4 )
+*           On entry ISEED specifies the seed of the random number
+*           generator. 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 DLATM7
+*           to continue the same random number sequence.
+*           Changed on exit.
+*
+*  D      - DOUBLE PRECISION array, dimension ( MIN( M , N ) )
+*           Array to be computed according to MODE, COND and IRSIGN.
+*           May be changed on exit if MODE is nonzero.
+*
+*  N      - INTEGER
+*           Number of entries of D. Not modified.
+*
+*  RANK   - INTEGER
+*           The rank of matrix to be generated for modes 1,2,3 only.
+*           D( RANK+1:N ) = 0.
+*           Not modified.
+*
+*  INFO   - INTEGER
+*            0  => normal termination
+*           -1  => if MODE not in range -6 to 6
+*           -2  => if MODE neither -6, 0 nor 6, and
+*                  IRSIGN neither 0 nor 1
+*           -3  => if MODE neither -6, 0 nor 6 and COND less than 1
+*           -4  => if MODE equals 6 or -6 and IDIST not in range 1 to 3
+*           -7  => if N negative
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D0 )
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   HALF
+      PARAMETER          ( HALF = 0.5D0 )
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION   ALPHA, TEMP
+      INTEGER            I
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLARAN
+      EXTERNAL           DLARAN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLARNV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DBLE, EXP, LOG
+*     ..
+*     .. Executable Statements ..
+*
+*     Decode and Test the input parameters. Initialize flags & seed.
+*
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Set INFO if an error
+*
+      IF( MODE.LT.-6 .OR. MODE.GT.6 ) THEN
+         INFO = -1
+      ELSE IF( ( MODE.NE.-6 .AND. MODE.NE.0 .AND. MODE.NE.6 ) .AND.
+     $         ( IRSIGN.NE.0 .AND. IRSIGN.NE.1 ) ) THEN
+         INFO = -2
+      ELSE IF( ( MODE.NE.-6 .AND. MODE.NE.0 .AND. MODE.NE.6 ) .AND.
+     $         COND.LT.ONE ) THEN
+         INFO = -3
+      ELSE IF( ( MODE.EQ.6 .OR. MODE.EQ.-6 ) .AND.
+     $         ( IDIST.LT.1 .OR. IDIST.GT.3 ) ) THEN
+         INFO = -4
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -7
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DLATM7', -INFO )
+         RETURN
+      END IF
+*
+*     Compute D according to COND and MODE
+*
+      IF( MODE.NE.0 ) THEN
+         GO TO ( 100, 130, 160, 190, 210, 230 )ABS( MODE )
+*
+*        One large D value:
+*
+  100    CONTINUE
+         DO 110 I = 2, RANK
+            D( I ) = ONE / COND
+  110    CONTINUE
+         DO 120 I = RANK + 1, N
+            D( I ) = ZERO
+  120    CONTINUE
+         D( 1 ) = ONE
+         GO TO 240
+*
+*        One small D value:
+*
+  130    CONTINUE
+         DO 140 I = 1, RANK - 1
+            D( I ) = ONE
+  140    CONTINUE
+         DO 150 I = RANK + 1, N
+            D( I ) = ZERO
+  150    CONTINUE
+         D( RANK ) = ONE / COND
+         GO TO 240
+*
+*        Exponentially distributed D values:
+*
+  160    CONTINUE
+         D( 1 ) = ONE
+         IF( N.GT.1 ) THEN
+            ALPHA = COND**( -ONE / DBLE( RANK-1 ) )
+            DO 170 I = 2, RANK
+               D( I ) = ALPHA**( I-1 )
+  170       CONTINUE
+            DO 180 I = RANK + 1, N
+               D( I ) = ZERO
+  180       CONTINUE
+         END IF
+         GO TO 240
+*
+*        Arithmetically distributed D values:
+*
+  190    CONTINUE
+         D( 1 ) = ONE
+         IF( N.GT.1 ) THEN
+            TEMP = ONE / COND
+            ALPHA = ( ONE-TEMP ) / DBLE( N-1 )
+            DO 200 I = 2, N
+               D( I ) = DBLE( N-I )*ALPHA + TEMP
+  200       CONTINUE
+         END IF
+         GO TO 240
+*
+*        Randomly distributed D values on ( 1/COND , 1):
+*
+  210    CONTINUE
+         ALPHA = LOG( ONE / COND )
+         DO 220 I = 1, N
+            D( I ) = EXP( ALPHA*DLARAN( ISEED ) )
+  220    CONTINUE
+         GO TO 240
+*
+*        Randomly distributed D values from IDIST
+*
+  230    CONTINUE
+         CALL DLARNV( IDIST, ISEED, N, D )
+*
+  240    CONTINUE
+*
+*        If MODE neither -6 nor 0 nor 6, and IRSIGN = 1, assign
+*        random signs to D
+*
+         IF( ( MODE.NE.-6 .AND. MODE.NE.0 .AND. MODE.NE.6 ) .AND.
+     $       IRSIGN.EQ.1 ) THEN
+            DO 250 I = 1, N
+               TEMP = DLARAN( ISEED )
+               IF( TEMP.GT.HALF )
+     $            D( I ) = -D( I )
+  250       CONTINUE
+         END IF
+*
+*        Reverse if MODE < 0
+*
+         IF( MODE.LT.0 ) THEN
+            DO 260 I = 1, N / 2
+               TEMP = D( I )
+               D( I ) = D( N+1-I )
+               D( N+1-I ) = TEMP
+  260       CONTINUE
+         END IF
+*
+      END IF
+*
+      RETURN
+*
+*     End of DLATM7
+*
+      END
diff --git a/TESTING/MATGEN/dlatmt.f b/TESTING/MATGEN/dlatmt.f
new file mode 100644
index 00000000..2ff32554
--- /dev/null
+++ b/TESTING/MATGEN/dlatmt.f
@@ -0,0 +1,1047 @@
+      SUBROUTINE DLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
+     $                   RANK, KL, KU, PACK, A, LDA, WORK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Craig Lucas, University of Manchester / NAG Ltd.
+*     October, 2008
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   COND, DMAX
+      INTEGER            INFO, KL, KU, LDA, M, MODE, N, RANK
+      CHARACTER          DIST, PACK, SYM
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), D( * ), WORK( * )
+      INTEGER            ISEED( 4 )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*     DLATMT generates random matrices with specified singular values
+*     (or symmetric/hermitian with specified eigenvalues)
+*     for testing LAPACK programs.
+*
+*     DLATMT 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.
+*
+*  Arguments
+*  =========
+*
+*  M      - INTEGER
+*           The number of rows of A. Not modified.
+*
+*  N      - INTEGER
+*           The number of columns of A. 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 DLATMT
+*           to continue the same random number sequence.
+*           Changed on exit.
+*
+*  SYM    - 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.
+*
+*  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:RANK)=1.0/COND
+*           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
+*           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
+*
+*           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
+*           MODE = 5 sets D to random numbers in the range
+*                    ( 1/COND , 1 ) such that their logarithms
+*                    are uniformly distributed.
+*           MODE = 6 set D to random numbers from same distribution
+*                    as the rest of the matrix.
+*           MODE < 0 has the same meaning as ABS(MODE), except that
+*              the order of the elements of D is reversed.
+*           Thus if MODE is positive, D has entries ranging from
+*              1 to 1/COND, if negative, from 1/COND to 1,
+*           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
+*              the elements of D will also be multiplied by a random
+*              sign (i.e., +1 or -1.)
+*           Not modified.
+*
+*  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.
+*
+*  RANK   - INTEGER
+*           The rank of matrix to be generated for modes 1,2,3 only.
+*           D( RANK+1:N ) = 0.
+*           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.
+*           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.
+*           Not modified.
+*
+*  PACK   - 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 DLATMT differ only in the PACK parameter,
+*           they will generate mathematically equivalent matrices.
+*           Not modified.
+*
+*  A      - 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.
+*
+*  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   - DOUBLE PRECISION 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='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 DLATM7
+*            2  => Cannot scale to DMAX (max. sing. value is 0)
+*            3  => Error return from DLAGGE or DLAGSY
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D0 )
+      DOUBLE PRECISION   TWOPI
+      PARAMETER          ( TWOPI = 6.2831853071795864769252867663D+0 )
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION   ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
+      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
+      LOGICAL            GIVENS, ILEXTR, ILTEMP, TOPDWN
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLARND
+      LOGICAL            LSAME
+      EXTERNAL           DLARND, LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLATM7, DCOPY, DLAGGE, DLAGSY, DLAROT,
+     $                   DLARTG, DLASET, DSCAL, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, COS, DBLE, 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
+      ELSE IF( LSAME( SYM, 'P' ) ) THEN
+         ISYM = 2
+         IRSIGN = 0
+      ELSE IF( LSAME( SYM, 'S' ) ) THEN
+         ISYM = 2
+         IRSIGN = 1
+      ELSE IF( LSAME( SYM, 'H' ) ) THEN
+         ISYM = 2
+         IRSIGN = 1
+      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 )
+*
+      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( 'DLATMT', -INFO )
+         RETURN
+      END IF
+*
+*     Initialize random number generator
+*
+      DO 100 I = 1, 4
+         ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
+  100 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 DLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, RANK,
+     $             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( RANK ) ) ) 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 110 I = 2, RANK
+            TEMP = MAX( TEMP, ABS( D( I ) ) )
+  110    CONTINUE
+*
+         IF( TEMP.GT.ZERO ) THEN
+            ALPHA = DMAX / TEMP
+         ELSE
+            INFO = 2
+            RETURN
+         END IF
+*
+         CALL DSCAL( RANK, ALPHA, D, 1 )
+*
+      END IF
+*
+*     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, SY, GB, 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
+      CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
+*
+*     Diagonal Matrix -- We are done, unless it
+*     is to be stored SP/PP/TP (PACK='R' or 'C')
+*
+      IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN
+         CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
+         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
+*
+            CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
+*
+            IF( TOPDWN ) THEN
+               JKL = 0
+               DO 140 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 130 JR = 1, MIN( M+JKU, N ) + JKL - 1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*DLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     ICOL = MAX( 1, JR-JKL )
+                     IF( JR.LT.M ) THEN
+                        IL = MIN( N, JR+JKU ) + 1 - ICOL
+                        CALL DLAROT( .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 120 JCH = JR - JKL, 1, -JKL - JKU
+                        IF( IR.LT.M ) THEN
+                           CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), EXTRA, C, S, DUMMY )
+                        END IF
+                        IROW = MAX( 1, JCH-JKU )
+                        IL = IR + 2 - IROW
+                        TEMP = ZERO
+                        ILTEMP = JCH.GT.JKU
+                        CALL DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S,
+     $                               A( IROW-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, TEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL DLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), TEMP, C, S, DUMMY )
+                           ICOL = MAX( 1, JCH-JKU-JKL )
+                           IL = IC + 2 - ICOL
+                           EXTRA = ZERO
+                           CALL DLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
+     $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
+     $                                  IOFFST, ICOL ), ILDA, EXTRA,
+     $                                  TEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  120                CONTINUE
+  130             CONTINUE
+  140          CONTINUE
+*
+               JKU = UUB
+               DO 170 JKL = 1, LLB
+*
+*                 Transform from bandwidth JKL-1, JKU to JKL, JKU
+*
+                  DO 160 JC = 1, MIN( N+JKL, M ) + JKU - 1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*DLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     IROW = MAX( 1, JC-JKU )
+                     IF( JC.LT.N ) THEN
+                        IL = MIN( M, JC+JKL ) + 1 - IROW
+                        CALL DLAROT( .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 150 JCH = JC - JKU, 1, -JKL - JKU
+                        IF( IC.LT.N ) THEN
+                           CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), EXTRA, C, S, DUMMY )
+                        END IF
+                        ICOL = MAX( 1, JCH-JKL )
+                        IL = IC + 2 - ICOL
+                        TEMP = ZERO
+                        ILTEMP = JCH.GT.JKL
+                        CALL DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S,
+     $                               A( IR-ISKEW*ICOL+IOFFST, ICOL ),
+     $                               ILDA, TEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL DLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
+     $                                  ICOL+1 ), TEMP, C, S, DUMMY )
+                           IROW = MAX( 1, JCH-JKL-JKU )
+                           IL = IR + 2 - IROW
+                           EXTRA = ZERO
+                           CALL DLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
+     $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
+     $                                  IOFFST, ICOL ), ILDA, EXTRA,
+     $                                  TEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  150                CONTINUE
+  160             CONTINUE
+  170          CONTINUE
+*
+            ELSE
+*
+*              Bottom-Up -- Start at the bottom right.
+*
+               JKL = 0
+               DO 200 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 190 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*DLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     IROW = MAX( 1, JC-JKU+1 )
+                     IF( JC.GT.0 ) THEN
+                        IL = MIN( M, JC+JKL+1 ) + 1 - IROW
+                        CALL DLAROT( .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 180 JCH = JC + JKL, IENDCH, JKL + JKU
+                        ILEXTR = IC.GT.0
+                        IF( ILEXTR ) THEN
+                           CALL DLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                                  EXTRA, C, S, DUMMY )
+                        END IF
+                        IC = MAX( 1, IC )
+                        ICOL = MIN( N-1, JCH+JKU )
+                        ILTEMP = JCH + JKU.LT.N
+                        TEMP = ZERO
+                        CALL DLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
+     $                               C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, EXTRA, TEMP )
+                        IF( ILTEMP ) THEN
+                           CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), TEMP, C, S, DUMMY )
+                           IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+                           EXTRA = ZERO
+                           CALL DLAROT( .FALSE., .TRUE.,
+     $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+     $                                  A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), ILDA, TEMP, EXTRA )
+                           IC = ICOL
+                        END IF
+  180                CONTINUE
+  190             CONTINUE
+  200          CONTINUE
+*
+               JKU = UUB
+               DO 230 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 220 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*DLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     ICOL = MAX( 1, JR-JKL+1 )
+                     IF( JR.GT.0 ) THEN
+                        IL = MIN( N, JR+JKU+1 ) + 1 - ICOL
+                        CALL DLAROT( .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 210 JCH = JR + JKU, IENDCH, JKL + JKU
+                        ILEXTR = IR.GT.0
+                        IF( ILEXTR ) THEN
+                           CALL DLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ),
+     $                                  EXTRA, C, S, DUMMY )
+                        END IF
+                        IR = MAX( 1, IR )
+                        IROW = MIN( M-1, JCH+JKL )
+                        ILTEMP = JCH + JKL.LT.M
+                        TEMP = ZERO
+                        CALL DLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
+     $                               C, S, A( IR-ISKEW*JCH+IOFFST,
+     $                               JCH ), ILDA, EXTRA, TEMP )
+                        IF( ILTEMP ) THEN
+                           CALL DLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  TEMP, C, S, DUMMY )
+                           IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+                           EXTRA = ZERO
+                           CALL DLAROT( .TRUE., .TRUE.,
+     $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+     $                                  A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  ILDA, TEMP, EXTRA )
+                           IR = IROW
+                        END IF
+  210                CONTINUE
+  220             CONTINUE
+  230          CONTINUE
+            END IF
+*
+         ELSE
+*
+*           Symmetric -- 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
+               CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
+*
+               DO 260 K = 1, UUB
+                  DO 250 JC = 1, N - 1
+                     IROW = MAX( 1, JC-K )
+                     IL = MIN( JC+1, K+2 )
+                     EXTRA = ZERO
+                     TEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 )
+                     ANGLE = TWOPI*DLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     CALL DLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
+     $                            A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
+     $                            EXTRA, TEMP )
+                     CALL DLAROT( .TRUE., .TRUE., .FALSE.,
+     $                            MIN( K, N-JC )+1, C, S,
+     $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+     $                            TEMP, DUMMY )
+*
+*                    Chase EXTRA back up the matrix
+*
+                     ICOL = JC
+                     DO 240 JCH = JC - K, 1, -K
+                        CALL DLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
+     $                               ICOL+1 ), EXTRA, C, S, DUMMY )
+                        TEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
+                        CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, -S,
+     $                               A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+     $                               ILDA, TEMP, EXTRA )
+                        IROW = MAX( 1, JCH-K )
+                        IL = MIN( JCH+1, K+2 )
+                        EXTRA = ZERO
+                        CALL DLAROT( .FALSE., JCH.GT.K, .TRUE., IL, C,
+     $                               -S, A( IROW-ISKEW*JCH+IOFFG, JCH ),
+     $                               ILDA, EXTRA, TEMP )
+                        ICOL = JCH
+  240                CONTINUE
+  250             CONTINUE
+  260          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 280 JC = 1, N
+                     IROW = IOFFST - ISKEW*JC
+                     DO 270 JR = JC, MIN( N, JC+UUB )
+                        A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  270                CONTINUE
+  280             CONTINUE
+                  IF( IPACK.EQ.5 ) THEN
+                     DO 300 JC = N - UUB + 1, N
+                        DO 290 JR = N + 2 - JC, UUB + 1
+                           A( JR, JC ) = ZERO
+  290                   CONTINUE
+  300                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
+               CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
+*
+               DO 330 K = 1, UUB
+                  DO 320 JC = N - 1, 1, -1
+                     IL = MIN( N+1-JC, K+2 )
+                     EXTRA = ZERO
+                     TEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
+                     ANGLE = TWOPI*DLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = -SIN( ANGLE )
+                     CALL DLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S,
+     $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+     $                            TEMP, EXTRA )
+                     ICOL = MAX( 1, JC-K+1 )
+                     CALL DLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, C,
+     $                            S, A( JC-ISKEW*ICOL+IOFFG, ICOL ),
+     $                            ILDA, DUMMY, TEMP )
+*
+*                    Chase EXTRA back down the matrix
+*
+                     ICOL = JC
+                     DO 310 JCH = JC + K, N - 1, K
+                        CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+     $                               EXTRA, C, S, DUMMY )
+                        TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
+                        CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
+     $                               A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+     $                               ILDA, EXTRA, TEMP )
+                        IL = MIN( N+1-JCH, K+2 )
+                        EXTRA = ZERO
+                        CALL DLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, C,
+     $                               S, A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+     $                               ILDA, TEMP, EXTRA )
+                        ICOL = JCH
+  310                CONTINUE
+  320             CONTINUE
+  330          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 350 JC = N, 1, -1
+                     IROW = IOFFST - ISKEW*JC
+                     DO 340 JR = JC, MAX( 1, JC-UUB ), -1
+                        A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  340                CONTINUE
+  350             CONTINUE
+                  IF( IPACK.EQ.6 ) THEN
+                     DO 370 JC = 1, UUB
+                        DO 360 JR = 1, UUB + 1 - JC
+                           A( JR, JC ) = ZERO
+  360                   CONTINUE
+  370                CONTINUE
+                  END IF
+                  IF( IPACKG.EQ.5 ) THEN
+                     IPACKG = IPACK
+                  ELSE
+                     IPACKG = 0
+                  END IF
+               END IF
+            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 DLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
+     $                   IINFO )
+         ELSE
+*
+*           Symmetric -- A = U D U'
+*
+            CALL DLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
+*
+         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 390 J = 1, M
+               DO 380 I = J + 1, M
+                  A( I, J ) = ZERO
+  380          CONTINUE
+  390       CONTINUE
+*
+         ELSE IF( IPACK.EQ.2 ) THEN
+*
+*           'L' -- Lower triangular, not packed
+*
+            DO 410 J = 2, M
+               DO 400 I = 1, J - 1
+                  A( I, J ) = ZERO
+  400          CONTINUE
+  410       CONTINUE
+*
+         ELSE IF( IPACK.EQ.3 ) THEN
+*
+*           'C' -- Upper triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 430 J = 1, M
+               DO 420 I = 1, J
+                  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.EQ.4 ) THEN
+*
+*           'R' -- Lower triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 450 J = 1, M
+               DO 440 I = J, M
+                  IROW = IROW + 1
+                  IF( IROW.GT.LDA ) THEN
+                     IROW = 1
+                     ICOL = ICOL + 1
+                  END IF
+                  A( IROW, ICOL ) = A( I, J )
+  440          CONTINUE
+  450       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 470 J = 1, UUB
+               DO 460 I = MIN( J+LLB, M ), 1, -1
+                  A( I-J+UUB+1, J ) = A( I, J )
+  460          CONTINUE
+  470       CONTINUE
+*
+            DO 490 J = UUB + 2, N
+               DO 480 I = J - UUB, MIN( J+LLB, M )
+                  A( I-J+UUB+1, J ) = A( I, J )
+  480          CONTINUE
+  490       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 510 JC = ICOL, M
+               DO 500 JR = IROW + 1, LDA
+                  A( JR, JC ) = ZERO
+  500          CONTINUE
+               IROW = 0
+  510       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 540 JC = 1, N
+               DO 520 JR = 1, UUB + 1 - JC
+                  A( JR, JC ) = ZERO
+  520          CONTINUE
+               DO 530 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
+                  A( JR, JC ) = ZERO
+  530          CONTINUE
+  540       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of DLATMT
+*
+      END
diff --git a/TESTING/MATGEN/slahilb.f b/TESTING/MATGEN/slahilb.f
new file mode 100644
index 00000000..afdc4201
--- /dev/null
+++ b/TESTING/MATGEN/slahilb.f
@@ -0,0 +1,166 @@
+      SUBROUTINE SLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK, INFO)
+!
+!  -- LAPACK auxiliary test routine (version 3.0) --
+!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+!     Courant Institute, Argonne National Lab, and Rice University
+!     28 August, 2006
+!
+!     David Vu <dtv@cs.berkeley.edu>      
+!     Yozo Hida <yozo@cs.berkeley.edu>      
+!     Jason Riedy <ejr@cs.berkeley.edu>
+!     D. Halligan <dhalligan@berkeley.edu>
+!
+      IMPLICIT NONE
+!     .. Scalar Arguments ..
+      INTEGER N, NRHS, LDA, LDX, LDB, INFO
+!     .. Array Arguments ..
+      REAL A(LDA, N), X(LDX, NRHS), B(LDB, NRHS), WORK(N)
+!     ..
+!
+!  Purpose
+!  =======
+!
+!  SLAHILB generates an N by N scaled Hilbert matrix in A along with
+!  NRHS right-hand sides in B and solutions in X such that A*X=B.
+!
+!  The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
+!  entries are integers.  The right-hand sides are the first NRHS 
+!  columns of M * the identity matrix, and the solutions are the 
+!  first NRHS columns of the inverse Hilbert matrix.
+!
+!  The condition number of the Hilbert matrix grows exponentially with
+!  its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
+!  Hilbert matrices beyond a relatively small dimension cannot be
+!  generated exactly without extra precision.  Precision is exhausted
+!  when the largest entry in the inverse Hilbert matrix is greater than
+!  2 to the power of the number of bits in the fraction of the data type
+!  used plus one, which is 24 for single precision.  
+!
+!  In single, the generated solution is exact for N <= 6 and has
+!  small componentwise error for 7 <= N <= 11.
+!
+!  Arguments
+!  =========
+!
+!  N       (input) INTEGER
+!          The dimension of the matrix A.
+!      
+!  NRHS    (input) NRHS
+!          The requested number of right-hand sides.
+!
+!  A       (output) REAL array, dimension (LDA, N)
+!          The generated scaled Hilbert matrix.
+!
+!  LDA     (input) INTEGER
+!          The leading dimension of the array A.  LDA >= N.
+!
+!  X       (output) REAL array, dimension (LDX, NRHS)
+!          The generated exact solutions.  Currently, the first NRHS
+!          columns of the inverse Hilbert matrix.
+!
+!  LDX     (input) INTEGER
+!          The leading dimension of the array X.  LDX >= N.
+!
+!  B       (output) REAL array, dimension (LDB, NRHS)
+!          The generated right-hand sides.  Currently, the first NRHS
+!          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
+!
+!  LDB     (input) INTEGER
+!          The leading dimension of the array B.  LDB >= N.
+!
+!  WORK    (workspace) REAL array, dimension (N)
+!
+!
+!  INFO    (output) INTEGER
+!          = 0: successful exit
+!          = 1: N is too large; the data is still generated but may not
+!               be not exact.
+!          < 0: if INFO = -i, the i-th argument had an illegal value
+!
+!  =====================================================================
+
+!     .. Local Scalars ..
+      INTEGER TM, TI, R
+      INTEGER M
+      INTEGER I, J
+
+!     .. Parameters ..
+!     NMAX_EXACT   the largest dimension where the generated data is
+!                  exact.
+!     NMAX_APPROX  the largest dimension where the generated data has
+!                  a small componentwise relative error.
+      INTEGER NMAX_EXACT, NMAX_APPROX
+      PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11)
+
+!     ..
+!     .. External Functions
+      EXTERNAL SLASET
+      INTRINSIC REAL
+!     ..
+!     .. Executable Statements ..
+!
+!     Test the input arguments
+!
+      INFO = 0
+      IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
+         INFO = -1
+      ELSE IF (NRHS .LT. 0) THEN
+         INFO = -2
+      ELSE IF (LDA .LT. N) THEN
+         INFO = -4
+      ELSE IF (LDX .LT. N) THEN
+         INFO = -6
+      ELSE IF (LDB .LT. N) THEN
+         INFO = -8
+      END IF
+      IF (INFO .LT. 0) THEN
+         CALL XERBLA('SLAHILB', -INFO)
+         RETURN
+      END IF
+      IF (N .GT. NMAX_EXACT) THEN
+         INFO = 1
+      END IF
+
+!     Compute M = the LCM of the integers [1, 2*N-1].  The largest
+!     reasonable N is small enough that integers suffice (up to N = 11).
+      M = 1
+      DO I = 2, (2*N-1)
+         TM = M
+         TI = I
+         R = MOD(TM, TI)
+         DO WHILE (R .NE. 0)
+            TM = TI
+            TI = R
+            R = MOD(TM, TI)
+         END DO
+         M = (M / TI) * I
+      END DO
+
+!     Generate the scaled Hilbert matrix in A
+      DO J = 1, N
+         DO I = 1, N
+            A(I, J) = REAL(M) / (I + J - 1)
+         END DO
+      END DO
+
+!     Generate matrix B as simply the first NRHS columns of M * the
+!     identity.
+      CALL SLASET('Full', N, NRHS, 0.0, REAL(M), B, LDB)
+
+!     Generate the true solutions in X.  Because B = the first NRHS
+!     columns of M*I, the true solutions are just the first NRHS columns
+!     of the inverse Hilbert matrix.
+      WORK(1) = N
+      DO J = 2, N
+         WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
+     $        * (N +J -1)
+      END DO
+      
+      DO J = 1, NRHS
+         DO I = 1, N
+            X(I, J) = (WORK(I)*WORK(J)) / (I + J - 1)
+         END DO
+      END DO
+
+      END
+
diff --git a/TESTING/MATGEN/slarot.f b/TESTING/MATGEN/slarot.f
index 694e9bb8..466cb29f 100644
--- a/TESTING/MATGEN/slarot.f
+++ b/TESTING/MATGEN/slarot.f
@@ -19,7 +19,6 @@
 *
 *     SLAROT applies a (Givens) rotation to two adjacent rows or
 *     columns, where one element of the first and/or last column/row
-*     November 2006
 *     for use on matrices stored in some format other than GE, so
 *     that elements of the matrix may be used or modified for which
 *     no array element is provided.
diff --git a/TESTING/MATGEN/slatm7.f b/TESTING/MATGEN/slatm7.f
new file mode 100644
index 00000000..8fa977f6
--- /dev/null
+++ b/TESTING/MATGEN/slatm7.f
@@ -0,0 +1,255 @@
+      SUBROUTINE SLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, N,
+     $                   RANK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Craig Lucas, University of Manchester / NAG Ltd.
+*     October, 2008
+*
+*     .. Scalar Arguments ..
+      REAL               COND
+      INTEGER            IDIST, INFO, IRSIGN, MODE, N, RANK
+*     ..
+*     .. Array Arguments ..
+      REAL               D( * )
+      INTEGER            ISEED( 4 )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*     SLATM7 computes the entries of D as specified by MODE
+*     COND and IRSIGN. IDIST and ISEED determine the generation
+*     of random numbers. SLATM7 is called by SLATMT to generate
+*     random test matrices.
+*
+*  Arguments
+*  =========
+*
+*  MODE   - INTEGER
+*           On entry describes how D is to be computed:
+*           MODE = 0 means do not change D.
+*
+*           MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
+*           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
+*           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1)) I=1:RANK
+*
+*           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,
+*           Not modified.
+*
+*  COND   - REAL
+*           On entry, used as described under MODE above.
+*           If used, it must be >= 1. Not modified.
+*
+*  IRSIGN - INTEGER
+*           On entry, if MODE neither -6, 0 nor 6, determines sign of
+*           entries of D
+*           0 => leave entries of D unchanged
+*           1 => multiply each entry of D by 1 or -1 with probability .5
+*
+*  IDIST  - CHARACTER*1
+*           On entry, IDIST specifies the type of distribution to be
+*           used to generate a random matrix .
+*           1 => UNIFORM( 0, 1 )
+*           2 => UNIFORM( -1, 1 )
+*           3 => NORMAL( 0, 1 )
+*           Not modified.
+*
+*  ISEED  - INTEGER array, dimension ( 4 )
+*           On entry ISEED specifies the seed of the random number
+*           generator. 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 SLATM7
+*           to continue the same random number sequence.
+*           Changed on exit.
+*
+*  D      - REAL array, dimension ( MIN( M , N ) )
+*           Array to be computed according to MODE, COND and IRSIGN.
+*           May be changed on exit if MODE is nonzero.
+*
+*  N      - INTEGER
+*           Number of entries of D. Not modified.
+*
+*  RANK   - INTEGER
+*           The rank of matrix to be generated for modes 1,2,3 only.
+*           D( RANK+1:N ) = 0.
+*           Not modified.
+*
+*  INFO   - INTEGER
+*            0  => normal termination
+*           -1  => if MODE not in range -6 to 6
+*           -2  => if MODE neither -6, 0 nor 6, and
+*                  IRSIGN neither 0 nor 1
+*           -3  => if MODE neither -6, 0 nor 6 and COND less than 1
+*           -4  => if MODE equals 6 or -6 and IDIST not in range 1 to 3
+*           -7  => if N negative
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE
+      PARAMETER          ( ONE = 1.0E0 )
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E0 )
+      REAL               HALF
+      PARAMETER          ( HALF = 0.5E0 )
+*     ..
+*     .. Local Scalars ..
+      REAL               ALPHA, TEMP
+      INTEGER            I
+*     ..
+*     .. External Functions ..
+      REAL               SLARAN
+      EXTERNAL           SLARAN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLARNV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, EXP, LOG, REAL
+*     ..
+*     .. Executable Statements ..
+*
+*     Decode and Test the input parameters. Initialize flags & seed.
+*
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Set INFO if an error
+*
+      IF( MODE.LT.-6 .OR. MODE.GT.6 ) THEN
+         INFO = -1
+      ELSE IF( ( MODE.NE.-6 .AND. MODE.NE.0 .AND. MODE.NE.6 ) .AND.
+     $         ( IRSIGN.NE.0 .AND. IRSIGN.NE.1 ) ) THEN
+         INFO = -2
+      ELSE IF( ( MODE.NE.-6 .AND. MODE.NE.0 .AND. MODE.NE.6 ) .AND.
+     $         COND.LT.ONE ) THEN
+         INFO = -3
+      ELSE IF( ( MODE.EQ.6 .OR. MODE.EQ.-6 ) .AND.
+     $         ( IDIST.LT.1 .OR. IDIST.GT.3 ) ) THEN
+         INFO = -4
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -7
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SLATM7', -INFO )
+         RETURN
+      END IF
+*
+*     Compute D according to COND and MODE
+*
+      IF( MODE.NE.0 ) THEN
+         GO TO ( 100, 130, 160, 190, 210, 230 )ABS( MODE )
+*
+*        One large D value:
+*
+  100    CONTINUE
+         DO 110 I = 2, RANK
+            D( I ) = ONE / COND
+  110    CONTINUE
+         DO 120 I = RANK + 1, N
+            D( I ) = ZERO
+  120    CONTINUE
+         D( 1 ) = ONE
+         GO TO 240
+*
+*        One small D value:
+*
+  130    CONTINUE
+         DO 140 I = 1, RANK - 1
+            D( I ) = ONE
+  140    CONTINUE
+         DO 150 I = RANK + 1, N
+            D( I ) = ZERO
+  150    CONTINUE
+         D( RANK ) = ONE / COND
+         GO TO 240
+*
+*        Exponentially distributed D values:
+*
+  160    CONTINUE
+         D( 1 ) = ONE
+         IF( N.GT.1 ) THEN
+            ALPHA = COND**( -ONE / REAL( RANK-1 ) )
+            DO 170 I = 2, RANK
+               D( I ) = ALPHA**( I-1 )
+  170       CONTINUE
+            DO 180 I = RANK + 1, N
+               D( I ) = ZERO
+  180       CONTINUE
+         END IF
+         GO TO 240
+*
+*        Arithmetically distributed D values:
+*
+  190    CONTINUE
+         D( 1 ) = ONE
+         IF( N.GT.1 ) THEN
+            TEMP = ONE / COND
+            ALPHA = ( ONE-TEMP ) / REAL( N-1 )
+            DO 200 I = 2, N
+               D( I ) = REAL( N-I )*ALPHA + TEMP
+  200       CONTINUE
+         END IF
+         GO TO 240
+*
+*        Randomly distributed D values on ( 1/COND , 1):
+*
+  210    CONTINUE
+         ALPHA = LOG( ONE / COND )
+         DO 220 I = 1, N
+            D( I ) = EXP( ALPHA*SLARAN( ISEED ) )
+  220    CONTINUE
+         GO TO 240
+*
+*        Randomly distributed D values from IDIST
+*
+  230    CONTINUE
+         CALL SLARNV( IDIST, ISEED, N, D )
+*
+  240    CONTINUE
+*
+*        If MODE neither -6 nor 0 nor 6, and IRSIGN = 1, assign
+*        random signs to D
+*
+         IF( ( MODE.NE.-6 .AND. MODE.NE.0 .AND. MODE.NE.6 ) .AND.
+     $       IRSIGN.EQ.1 ) THEN
+            DO 250 I = 1, N
+               TEMP = SLARAN( ISEED )
+               IF( TEMP.GT.HALF )
+     $            D( I ) = -D( I )
+  250       CONTINUE
+         END IF
+*
+*        Reverse if MODE < 0
+*
+         IF( MODE.LT.0 ) THEN
+            DO 260 I = 1, N / 2
+               TEMP = D( I )
+               D( I ) = D( N+1-I )
+               D( N+1-I ) = TEMP
+  260       CONTINUE
+         END IF
+*
+      END IF
+*
+      RETURN
+*
+*     End of SLATM7
+*
+      END
diff --git a/TESTING/MATGEN/slatmt.f b/TESTING/MATGEN/slatmt.f
new file mode 100644
index 00000000..c0119cc5
--- /dev/null
+++ b/TESTING/MATGEN/slatmt.f
@@ -0,0 +1,1047 @@
+      SUBROUTINE SLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
+     $                   RANK, KL, KU, PACK, A, LDA, WORK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Craig Lucas, University of Manchester / NAG Ltd.
+*     October, 2008
+*
+*     .. Scalar Arguments ..
+      REAL               COND, DMAX
+      INTEGER            INFO, KL, KU, LDA, M, MODE, N, RANK
+      CHARACTER          DIST, PACK, SYM
+*     ..
+*     .. Array Arguments ..
+      REAL               A( LDA, * ), D( * ), WORK( * )
+      INTEGER            ISEED( 4 )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*     SLATMT generates random matrices with specified singular values
+*     (or symmetric/hermitian with specified eigenvalues)
+*     for testing LAPACK programs.
+*
+*     SLATMT operates by applying the following sequence of
+*     operations:
+*
+*       Set the diagonal to D, where D may be input or
+*          computed according to MODE, COND, DMAX, and SYM
+*          as described below.
+*
+*       Generate a matrix with the appropriate band structure, by one
+*          of two methods:
+*
+*       Method A:
+*           Generate a dense M x N matrix by multiplying D on the left
+*               and the right by random unitary matrices, then:
+*
+*           Reduce the bandwidth according to KL and KU, using
+*           Householder transformations.
+*
+*       Method B:
+*           Convert the bandwidth-0 (i.e., diagonal) matrix to a
+*               bandwidth-1 matrix using Givens rotations, "chasing"
+*               out-of-band elements back, much as in QR; then
+*               convert the bandwidth-1 to a bandwidth-2 matrix, etc.
+*               Note that for reasonably small bandwidths (relative to
+*               M and N) this requires less storage, as a dense matrix
+*               is not generated.  Also, for symmetric matrices, only
+*               one triangle is generated.
+*
+*       Method A is chosen if the bandwidth is a large fraction of the
+*           order of the matrix, and LDA is at least M (so a dense
+*           matrix can be stored.)  Method B is chosen if the bandwidth
+*           is small (< 1/2 N for symmetric, < .3 N+M for
+*           non-symmetric), or LDA is less than M and not less than the
+*           bandwidth.
+*
+*       Pack the matrix if desired. Options specified by PACK are:
+*          no packing
+*          zero out upper half (if symmetric)
+*          zero out lower half (if symmetric)
+*          store the upper half columnwise (if symmetric or upper
+*                triangular)
+*          store the lower half columnwise (if symmetric or lower
+*                triangular)
+*          store the lower triangle in banded format (if symmetric
+*                or lower triangular)
+*          store the upper triangle in banded format (if symmetric
+*                or upper triangular)
+*          store the entire matrix in banded format
+*       If Method B is chosen, and band format is specified, then the
+*          matrix will be generated in the band format, so no repacking
+*          will be necessary.
+*
+*  Arguments
+*  =========
+*
+*  M      - INTEGER
+*           The number of rows of A. Not modified.
+*
+*  N      - INTEGER
+*           The number of columns of A. 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 SLATMT
+*           to continue the same random number sequence.
+*           Changed on exit.
+*
+*  SYM    - 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.
+*
+*  D      - REAL array, dimension ( MIN( M , N ) )
+*           This array is used to specify the singular values or
+*           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
+*           assumed to contain the singular/eigenvalues, otherwise
+*           they will be computed according to MODE, COND, and DMAX,
+*           and placed in D.
+*           Modified if MODE is nonzero.
+*
+*  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:RANK)=1.0/COND
+*           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
+*           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
+*
+*           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
+*           MODE = 5 sets D to random numbers in the range
+*                    ( 1/COND , 1 ) such that their logarithms
+*                    are uniformly distributed.
+*           MODE = 6 set D to random numbers from same distribution
+*                    as the rest of the matrix.
+*           MODE < 0 has the same meaning as ABS(MODE), except that
+*              the order of the elements of D is reversed.
+*           Thus if MODE is positive, D has entries ranging from
+*              1 to 1/COND, if negative, from 1/COND to 1,
+*           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
+*              the elements of D will also be multiplied by a random
+*              sign (i.e., +1 or -1.)
+*           Not modified.
+*
+*  COND   - REAL
+*           On entry, this is used as described under MODE above.
+*           If used, it must be >= 1. Not modified.
+*
+*  DMAX   - REAL
+*           If MODE is neither -6, 0 nor 6, the contents of D, as
+*           computed according to MODE and COND, will be scaled by
+*           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
+*           singular value (which is to say the norm) will be abs(DMAX).
+*           Note that DMAX need not be positive: if DMAX is negative
+*           (or zero), D will be scaled by a negative number (or zero).
+*           Not modified.
+*
+*  RANK   - INTEGER
+*           The rank of matrix to be generated for modes 1,2,3 only.
+*           D( RANK+1:N ) = 0.
+*           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.
+*           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.
+*           Not modified.
+*
+*  PACK   - CHARACTER*1
+*           This specifies packing of matrix as follows:
+*           'N' => no packing
+*           'U' => zero out all subdiagonal entries (if symmetric)
+*           'L' => zero out all superdiagonal entries (if symmetric)
+*           'C' => store the upper triangle columnwise
+*                  (only if the matrix is symmetric or upper triangular)
+*           'R' => store the lower triangle columnwise
+*                  (only if the matrix is symmetric or lower triangular)
+*           'B' => store the lower triangle in band storage scheme
+*                  (only if matrix symmetric or lower triangular)
+*           'Q' => store the upper triangle in band storage scheme
+*                  (only if matrix symmetric or upper triangular)
+*           'Z' => store the entire matrix in band storage scheme
+*                      (pivoting can be provided for by using this
+*                      option to store A in the trailing rows of
+*                      the allocated storage)
+*
+*           Using these options, the various LAPACK packed and banded
+*           storage schemes can be obtained:
+*           GB               - use 'Z'
+*           PB, SB or TB     - use 'B' or 'Q'
+*           PP, SP or TP     - use 'C' or 'R'
+*
+*           If two calls to SLATMT differ only in the PACK parameter,
+*           they will generate mathematically equivalent matrices.
+*           Not modified.
+*
+*  A      - REAL array, dimension ( LDA, N )
+*           On exit A is the desired test matrix.  A is first generated
+*           in full (unpacked) form, and then packed, if so specified
+*           by PACK.  Thus, the first M elements of the first N
+*           columns will always be modified.  If PACK specifies a
+*           packed or banded storage scheme, all LDA elements of the
+*           first N columns will be modified; the elements of the
+*           array which do not correspond to elements of the generated
+*           matrix are set to zero.
+*           Modified.
+*
+*  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   - REAL 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='S' or 'H' and KU not equal to KL
+*           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
+*                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
+*                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
+*                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
+*                  N.
+*           -14 => LDA is less than M, or PACK='Z' and LDA is less than
+*                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
+*            1  => Error return from SLATM7
+*            2  => Cannot scale to DMAX (max. sing. value is 0)
+*            3  => Error return from SLAGGE or SLAGSY
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E0 )
+      REAL               ONE
+      PARAMETER          ( ONE = 1.0E0 )
+      REAL               TWOPI
+      PARAMETER          ( TWOPI = 6.2831853071795864769252867663E+0 )
+*     ..
+*     .. Local Scalars ..
+      REAL               ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
+      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
+      LOGICAL            GIVENS, ILEXTR, ILTEMP, TOPDWN
+*     ..
+*     .. External Functions ..
+      REAL               SLARND
+      LOGICAL            LSAME
+      EXTERNAL           SLARND, LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLATM7, SCOPY, SLAGGE, SLAGSY, SLAROT,
+     $                   SLARTG, SLASET, SSCAL, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, COS, MAX, MIN, MOD, REAL, 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
+      ELSE IF( LSAME( SYM, 'P' ) ) THEN
+         ISYM = 2
+         IRSIGN = 0
+      ELSE IF( LSAME( SYM, 'S' ) ) THEN
+         ISYM = 2
+         IRSIGN = 1
+      ELSE IF( LSAME( SYM, 'H' ) ) THEN
+         ISYM = 2
+         IRSIGN = 1
+      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 )
+*
+      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( REAL( LLB+UUB ).LT.0.3*REAL( 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( 'SLATMT', -INFO )
+         RETURN
+      END IF
+*
+*     Initialize random number generator
+*
+      DO 100 I = 1, 4
+         ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
+  100 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 SLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, RANK,
+     $             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( RANK ) ) ) 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 110 I = 2, RANK
+            TEMP = MAX( TEMP, ABS( D( I ) ) )
+  110    CONTINUE
+*
+         IF( TEMP.GT.ZERO ) THEN
+            ALPHA = DMAX / TEMP
+         ELSE
+            INFO = 2
+            RETURN
+         END IF
+*
+         CALL SSCAL( RANK, ALPHA, D, 1 )
+*
+      END IF
+*
+*     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, SY, GB, 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
+      CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
+*
+*     Diagonal Matrix -- We are done, unless it
+*     is to be stored SP/PP/TP (PACK='R' or 'C')
+*
+      IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN
+         CALL SCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
+         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
+*
+            CALL SCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
+*
+            IF( TOPDWN ) THEN
+               JKL = 0
+               DO 140 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 130 JR = 1, MIN( M+JKU, N ) + JKL - 1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     ICOL = MAX( 1, JR-JKL )
+                     IF( JR.LT.M ) THEN
+                        IL = MIN( N, JR+JKU ) + 1 - ICOL
+                        CALL SLAROT( .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 120 JCH = JR - JKL, 1, -JKL - JKU
+                        IF( IR.LT.M ) THEN
+                           CALL SLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), EXTRA, C, S, DUMMY )
+                        END IF
+                        IROW = MAX( 1, JCH-JKU )
+                        IL = IR + 2 - IROW
+                        TEMP = ZERO
+                        ILTEMP = JCH.GT.JKU
+                        CALL SLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S,
+     $                               A( IROW-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, TEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL SLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), TEMP, C, S, DUMMY )
+                           ICOL = MAX( 1, JCH-JKU-JKL )
+                           IL = IC + 2 - ICOL
+                           EXTRA = ZERO
+                           CALL SLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
+     $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
+     $                                  IOFFST, ICOL ), ILDA, EXTRA,
+     $                                  TEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  120                CONTINUE
+  130             CONTINUE
+  140          CONTINUE
+*
+               JKU = UUB
+               DO 170 JKL = 1, LLB
+*
+*                 Transform from bandwidth JKL-1, JKU to JKL, JKU
+*
+                  DO 160 JC = 1, MIN( N+JKL, M ) + JKU - 1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     IROW = MAX( 1, JC-JKU )
+                     IF( JC.LT.N ) THEN
+                        IL = MIN( M, JC+JKL ) + 1 - IROW
+                        CALL SLAROT( .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 150 JCH = JC - JKU, 1, -JKL - JKU
+                        IF( IC.LT.N ) THEN
+                           CALL SLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), EXTRA, C, S, DUMMY )
+                        END IF
+                        ICOL = MAX( 1, JCH-JKL )
+                        IL = IC + 2 - ICOL
+                        TEMP = ZERO
+                        ILTEMP = JCH.GT.JKL
+                        CALL SLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S,
+     $                               A( IR-ISKEW*ICOL+IOFFST, ICOL ),
+     $                               ILDA, TEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL SLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
+     $                                  ICOL+1 ), TEMP, C, S, DUMMY )
+                           IROW = MAX( 1, JCH-JKL-JKU )
+                           IL = IR + 2 - IROW
+                           EXTRA = ZERO
+                           CALL SLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
+     $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
+     $                                  IOFFST, ICOL ), ILDA, EXTRA,
+     $                                  TEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  150                CONTINUE
+  160             CONTINUE
+  170          CONTINUE
+*
+            ELSE
+*
+*              Bottom-Up -- Start at the bottom right.
+*
+               JKL = 0
+               DO 200 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 190 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     IROW = MAX( 1, JC-JKU+1 )
+                     IF( JC.GT.0 ) THEN
+                        IL = MIN( M, JC+JKL+1 ) + 1 - IROW
+                        CALL SLAROT( .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 180 JCH = JC + JKL, IENDCH, JKL + JKU
+                        ILEXTR = IC.GT.0
+                        IF( ILEXTR ) THEN
+                           CALL SLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                                  EXTRA, C, S, DUMMY )
+                        END IF
+                        IC = MAX( 1, IC )
+                        ICOL = MIN( N-1, JCH+JKU )
+                        ILTEMP = JCH + JKU.LT.N
+                        TEMP = ZERO
+                        CALL SLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
+     $                               C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, EXTRA, TEMP )
+                        IF( ILTEMP ) THEN
+                           CALL SLARTG( A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), TEMP, C, S, DUMMY )
+                           IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+                           EXTRA = ZERO
+                           CALL SLAROT( .FALSE., .TRUE.,
+     $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+     $                                  A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), ILDA, TEMP, EXTRA )
+                           IC = ICOL
+                        END IF
+  180                CONTINUE
+  190             CONTINUE
+  200          CONTINUE
+*
+               JKU = UUB
+               DO 230 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 220 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
+                     EXTRA = ZERO
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     ICOL = MAX( 1, JR-JKL+1 )
+                     IF( JR.GT.0 ) THEN
+                        IL = MIN( N, JR+JKU+1 ) + 1 - ICOL
+                        CALL SLAROT( .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 210 JCH = JR + JKU, IENDCH, JKL + JKU
+                        ILEXTR = IR.GT.0
+                        IF( ILEXTR ) THEN
+                           CALL SLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ),
+     $                                  EXTRA, C, S, DUMMY )
+                        END IF
+                        IR = MAX( 1, IR )
+                        IROW = MIN( M-1, JCH+JKL )
+                        ILTEMP = JCH + JKL.LT.M
+                        TEMP = ZERO
+                        CALL SLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
+     $                               C, S, A( IR-ISKEW*JCH+IOFFST,
+     $                               JCH ), ILDA, EXTRA, TEMP )
+                        IF( ILTEMP ) THEN
+                           CALL SLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  TEMP, C, S, DUMMY )
+                           IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
+                           EXTRA = ZERO
+                           CALL SLAROT( .TRUE., .TRUE.,
+     $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
+     $                                  A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  ILDA, TEMP, EXTRA )
+                           IR = IROW
+                        END IF
+  210                CONTINUE
+  220             CONTINUE
+  230          CONTINUE
+            END IF
+*
+         ELSE
+*
+*           Symmetric -- 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
+               CALL SCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
+*
+               DO 260 K = 1, UUB
+                  DO 250 JC = 1, N - 1
+                     IROW = MAX( 1, JC-K )
+                     IL = MIN( JC+1, K+2 )
+                     EXTRA = ZERO
+                     TEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 )
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = SIN( ANGLE )
+                     CALL SLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
+     $                            A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
+     $                            EXTRA, TEMP )
+                     CALL SLAROT( .TRUE., .TRUE., .FALSE.,
+     $                            MIN( K, N-JC )+1, C, S,
+     $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+     $                            TEMP, DUMMY )
+*
+*                    Chase EXTRA back up the matrix
+*
+                     ICOL = JC
+                     DO 240 JCH = JC - K, 1, -K
+                        CALL SLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
+     $                               ICOL+1 ), EXTRA, C, S, DUMMY )
+                        TEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
+                        CALL SLAROT( .TRUE., .TRUE., .TRUE., K+2, C, -S,
+     $                               A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+     $                               ILDA, TEMP, EXTRA )
+                        IROW = MAX( 1, JCH-K )
+                        IL = MIN( JCH+1, K+2 )
+                        EXTRA = ZERO
+                        CALL SLAROT( .FALSE., JCH.GT.K, .TRUE., IL, C,
+     $                               -S, A( IROW-ISKEW*JCH+IOFFG, JCH ),
+     $                               ILDA, EXTRA, TEMP )
+                        ICOL = JCH
+  240                CONTINUE
+  250             CONTINUE
+  260          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 280 JC = 1, N
+                     IROW = IOFFST - ISKEW*JC
+                     DO 270 JR = JC, MIN( N, JC+UUB )
+                        A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  270                CONTINUE
+  280             CONTINUE
+                  IF( IPACK.EQ.5 ) THEN
+                     DO 300 JC = N - UUB + 1, N
+                        DO 290 JR = N + 2 - JC, UUB + 1
+                           A( JR, JC ) = ZERO
+  290                   CONTINUE
+  300                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
+               CALL SCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
+*
+               DO 330 K = 1, UUB
+                  DO 320 JC = N - 1, 1, -1
+                     IL = MIN( N+1-JC, K+2 )
+                     EXTRA = ZERO
+                     TEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
+                     ANGLE = TWOPI*SLARND( 1, ISEED )
+                     C = COS( ANGLE )
+                     S = -SIN( ANGLE )
+                     CALL SLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S,
+     $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+     $                            TEMP, EXTRA )
+                     ICOL = MAX( 1, JC-K+1 )
+                     CALL SLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, C,
+     $                            S, A( JC-ISKEW*ICOL+IOFFG, ICOL ),
+     $                            ILDA, DUMMY, TEMP )
+*
+*                    Chase EXTRA back down the matrix
+*
+                     ICOL = JC
+                     DO 310 JCH = JC + K, N - 1, K
+                        CALL SLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+     $                               EXTRA, C, S, DUMMY )
+                        TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
+                        CALL SLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
+     $                               A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
+     $                               ILDA, EXTRA, TEMP )
+                        IL = MIN( N+1-JCH, K+2 )
+                        EXTRA = ZERO
+                        CALL SLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, C,
+     $                               S, A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+     $                               ILDA, TEMP, EXTRA )
+                        ICOL = JCH
+  310                CONTINUE
+  320             CONTINUE
+  330          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 350 JC = N, 1, -1
+                     IROW = IOFFST - ISKEW*JC
+                     DO 340 JR = JC, MAX( 1, JC-UUB ), -1
+                        A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  340                CONTINUE
+  350             CONTINUE
+                  IF( IPACK.EQ.6 ) THEN
+                     DO 370 JC = 1, UUB
+                        DO 360 JR = 1, UUB + 1 - JC
+                           A( JR, JC ) = ZERO
+  360                   CONTINUE
+  370                CONTINUE
+                  END IF
+                  IF( IPACKG.EQ.5 ) THEN
+                     IPACKG = IPACK
+                  ELSE
+                     IPACKG = 0
+                  END IF
+               END IF
+            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 SLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
+     $                   IINFO )
+         ELSE
+*
+*           Symmetric -- A = U D U'
+*
+            CALL SLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
+*
+         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 390 J = 1, M
+               DO 380 I = J + 1, M
+                  A( I, J ) = ZERO
+  380          CONTINUE
+  390       CONTINUE
+*
+         ELSE IF( IPACK.EQ.2 ) THEN
+*
+*           'L' -- Lower triangular, not packed
+*
+            DO 410 J = 2, M
+               DO 400 I = 1, J - 1
+                  A( I, J ) = ZERO
+  400          CONTINUE
+  410       CONTINUE
+*
+         ELSE IF( IPACK.EQ.3 ) THEN
+*
+*           'C' -- Upper triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 430 J = 1, M
+               DO 420 I = 1, J
+                  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.EQ.4 ) THEN
+*
+*           'R' -- Lower triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 450 J = 1, M
+               DO 440 I = J, M
+                  IROW = IROW + 1
+                  IF( IROW.GT.LDA ) THEN
+                     IROW = 1
+                     ICOL = ICOL + 1
+                  END IF
+                  A( IROW, ICOL ) = A( I, J )
+  440          CONTINUE
+  450       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 470 J = 1, UUB
+               DO 460 I = MIN( J+LLB, M ), 1, -1
+                  A( I-J+UUB+1, J ) = A( I, J )
+  460          CONTINUE
+  470       CONTINUE
+*
+            DO 490 J = UUB + 2, N
+               DO 480 I = J - UUB, MIN( J+LLB, M )
+                  A( I-J+UUB+1, J ) = A( I, J )
+  480          CONTINUE
+  490       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 510 JC = ICOL, M
+               DO 500 JR = IROW + 1, LDA
+                  A( JR, JC ) = ZERO
+  500          CONTINUE
+               IROW = 0
+  510       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 540 JC = 1, N
+               DO 520 JR = 1, UUB + 1 - JC
+                  A( JR, JC ) = ZERO
+  520          CONTINUE
+               DO 530 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
+                  A( JR, JC ) = ZERO
+  530          CONTINUE
+  540       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of SLATMT
+*
+      END
diff --git a/TESTING/MATGEN/zlahilb.f b/TESTING/MATGEN/zlahilb.f
new file mode 100644
index 00000000..9350684b
--- /dev/null
+++ b/TESTING/MATGEN/zlahilb.f
@@ -0,0 +1,206 @@
+      SUBROUTINE ZLAHILB(N, NRHS, A, LDA, X, LDX, B, LDB, WORK,
+     $     INFO, PATH)
+!
+!  -- LAPACK auxiliary test routine (version 3.0) --
+!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+!     Courant Institute, Argonne National Lab, and Rice University
+!     28 August, 2006
+!
+!     David Vu <dtv@cs.berkeley.edu>
+!     Yozo Hida <yozo@cs.berkeley.edu>
+!     Jason Riedy <ejr@cs.berkeley.edu>
+!     D. Halligan <dhalligan@berkeley.edu>
+!
+      IMPLICIT NONE
+!     .. Scalar Arguments ..
+      INTEGER N, NRHS, LDA, LDX, LDB, INFO
+!     .. Array Arguments ..
+      DOUBLE PRECISION WORK(N)
+      COMPLEX*16 A(LDA,N), X(LDX, NRHS), B(LDB, NRHS)
+      CHARACTER*3 PATH
+!     ..
+!
+!  Purpose
+!  =======
+!
+!  ZLAHILB generates an N by N scaled Hilbert matrix in A along with
+!  NRHS right-hand sides in B and solutions in X such that A*X=B.
+!
+!  The Hilbert matrix is scaled by M = LCM(1, 2, ..., 2*N-1) so that all
+!  entries are integers.  The right-hand sides are the first NRHS
+!  columns of M * the identity matrix, and the solutions are the
+!  first NRHS columns of the inverse Hilbert matrix.
+!
+!  The condition number of the Hilbert matrix grows exponentially with
+!  its size, roughly as O(e ** (3.5*N)).  Additionally, the inverse
+!  Hilbert matrices beyond a relatively small dimension cannot be
+!  generated exactly without extra precision.  Precision is exhausted
+!  when the largest entry in the inverse Hilbert matrix is greater than
+!  2 to the power of the number of bits in the fraction of the data type
+!  used plus one, which is 24 for single precision.
+!
+!  In single, the generated solution is exact for N <= 6 and has
+!  small componentwise error for 7 <= N <= 11.
+!
+!  Arguments
+!  =========
+!
+!  N       (input) INTEGER
+!          The dimension of the matrix A.
+!
+!  NRHS    (input) NRHS
+!          The requested number of right-hand sides.
+!
+!  A       (output) COMPLEX array, dimension (LDA, N)
+!          The generated scaled Hilbert matrix.
+!
+!  LDA     (input) INTEGER
+!          The leading dimension of the array A.  LDA >= N.
+!
+!  X       (output) COMPLEX array, dimension (LDX, NRHS)
+!          The generated exact solutions.  Currently, the first NRHS
+!          columns of the inverse Hilbert matrix.
+!
+!  LDX     (input) INTEGER
+!          The leading dimension of the array X.  LDX >= N.
+!
+!  B       (output) REAL array, dimension (LDB, NRHS)
+!          The generated right-hand sides.  Currently, the first NRHS
+!          columns of LCM(1, 2, ..., 2*N-1) * the identity matrix.
+!
+!  LDB     (input) INTEGER
+!          The leading dimension of the array B.  LDB >= N.
+!
+!  WORK    (workspace) REAL array, dimension (N)
+!
+!
+!  INFO    (output) INTEGER
+!          = 0: successful exit
+!          = 1: N is too large; the data is still generated but may not
+!               be not exact.
+!          < 0: if INFO = -i, the i-th argument had an illegal value
+!
+!  =====================================================================
+
+!     .. Local Scalars ..
+      INTEGER TM, TI, R
+      INTEGER M
+      INTEGER I, J
+      COMPLEX*16 TMP
+      CHARACTER*2 C2
+
+!     .. Parameters ..
+!     NMAX_EXACT   the largest dimension where the generated data is
+!                  exact.
+!     NMAX_APPROX  the largest dimension where the generated data has
+!                  a small componentwise relative error.
+!     ??? complex uses how many bits ???
+      INTEGER NMAX_EXACT, NMAX_APPROX, SIZE_D
+      PARAMETER (NMAX_EXACT = 6, NMAX_APPROX = 11, SIZE_D = 8)
+
+!     d's are generated from random permuation of those eight elements.
+      COMPLEX*16 d1(8), d2(8), invd1(8), invd2(8)
+      DATA D1 /(-1,0),(0,1),(-1,-1),(0,-1),(1,0),(-1,1),(1,1),(1,-1)/
+      DATA D2 /(-1,0),(0,-1),(-1,1),(0,1),(1,0),(-1,-1),(1,-1),(1,1)/
+      
+      DATA INVD1 /(-1,0),(0,-1),(-.5,.5),(0,1),(1,0),
+     $     (-.5,-.5),(.5,-.5),(.5,.5)/
+      DATA INVD2 /(-1,0),(0,1),(-.5,-.5),(0,-1),(1,0),
+     $     (-.5,.5),(.5,.5),(.5,-.5)/
+!     ..
+!     .. External Functions
+      EXTERNAL ZLASET, LSAMEN
+      INTRINSIC DBLE
+      LOGICAL LSAMEN
+!     ..
+!     .. Executable Statements ..
+      C2 = PATH( 2: 3 )
+!
+!     Test the input arguments
+!
+      INFO = 0
+      IF (N .LT. 0 .OR. N .GT. NMAX_APPROX) THEN
+         INFO = -1
+      ELSE IF (NRHS .LT. 0) THEN
+         INFO = -2
+      ELSE IF (LDA .LT. N) THEN
+         INFO = -4
+      ELSE IF (LDX .LT. N) THEN
+         INFO = -6
+      ELSE IF (LDB .LT. N) THEN
+         INFO = -8
+      END IF
+      IF (INFO .LT. 0) THEN
+         CALL XERBLA('ZLAHILB', -INFO)
+         RETURN
+      END IF
+      IF (N .GT. NMAX_EXACT) THEN
+         INFO = 1
+      END IF
+
+!     Compute M = the LCM of the integers [1, 2*N-1].  The largest
+!     reasonable N is small enough that integers suffice (up to N = 11).
+      M = 1
+      DO I = 2, (2*N-1)
+         TM = M
+         TI = I
+         R = MOD(TM, TI)
+         DO WHILE (R .NE. 0)
+            TM = TI
+            TI = R
+            R = MOD(TM, TI)
+         END DO
+         M = (M / TI) * I
+      END DO
+
+!     Generate the scaled Hilbert matrix in A
+!     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i*
+      IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
+         DO J = 1, N
+            DO I = 1, N
+               A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1))
+     $              * D1(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      ELSE
+         DO J = 1, N
+            DO I = 1, N
+               A(I, J) = D1(MOD(J,SIZE_D)+1) * (DBLE(M) / (I + J - 1))
+     $              * D2(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      END IF
+      
+!     Generate matrix B as simply the first NRHS columns of M * the
+!     identity.
+      TMP = DBLE(M)
+      CALL ZLASET('Full', N, NRHS, (0.0D+0,0.0D+0), TMP, B, LDB)
+
+!     Generate the true solutions in X.  Because B = the first NRHS
+!     columns of M*I, the true solutions are just the first NRHS columns
+!     of the inverse Hilbert matrix.
+      WORK(1) = N
+      DO J = 2, N
+         WORK(J) = (  ( (WORK(J-1)/(J-1)) * (J-1 - N) ) /(J-1)  )
+     $        * (N +J -1)
+      END DO
+
+!     If we are testing SY routines, take D1_i = D2_i, else, D1_i = D2_i*
+      IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
+         DO J = 1, NRHS
+            DO I = 1, N
+               X(I, J) = INVD1(MOD(J,SIZE_D)+1) *
+     $              ((WORK(I)*WORK(J)) / (I + J - 1))
+     $              * INVD1(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      ELSE
+         DO J = 1, NRHS
+            DO I = 1, N
+               X(I, J) = INVD2(MOD(J,SIZE_D)+1) *
+     $              ((WORK(I)*WORK(J)) / (I + J - 1))
+     $              * INVD1(MOD(I,SIZE_D)+1)
+            END DO
+         END DO
+      END IF
+      END
diff --git a/TESTING/MATGEN/zlarot.f b/TESTING/MATGEN/zlarot.f
index a5fec614..106a4d9f 100644
--- a/TESTING/MATGEN/zlarot.f
+++ b/TESTING/MATGEN/zlarot.f
@@ -19,7 +19,6 @@
 *
 *     ZLAROT applies a (Givens) rotation to two adjacent rows or
 *     columns, where one element of the first and/or last column/row
-*     November 2006
 *     for use on matrices stored in some format other than GE, so
 *     that elements of the matrix may be used or modified for which
 *     no array element is provided.
diff --git a/TESTING/MATGEN/zlatmt.f b/TESTING/MATGEN/zlatmt.f
new file mode 100644
index 00000000..69d6228b
--- /dev/null
+++ b/TESTING/MATGEN/zlatmt.f
@@ -0,0 +1,1172 @@
+      SUBROUTINE ZLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
+     $                   RANK, KL, KU, PACK, A, LDA, WORK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Craig Lucas, University of Manchester / NAG Ltd.
+*     October, 2008
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   COND, DMAX
+      INTEGER            INFO, KL, KU, LDA, M, MODE, N, RANK
+      CHARACTER          DIST, PACK, SYM
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16         A( LDA, * ), WORK( * )
+      DOUBLE PRECISION   D( * )
+      INTEGER            ISEED( 4 )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*     ZLATMT generates random matrices with specified singular values
+*     (or hermitian with specified eigenvalues)
+*     for testing LAPACK programs.
+*
+*     ZLATMT 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 ZLATMT
+*           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:RANK)=1.0/COND
+*           MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
+*           MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
+*           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
+*           MODE = 5 sets D to random numbers in the range
+*                    ( 1/COND , 1 ) such that their logarithms
+*                    are uniformly distributed.
+*           MODE = 6 set D to random numbers from same distribution
+*                    as the rest of the matrix.
+*           MODE < 0 has the same meaning as ABS(MODE), except that
+*              the order of the elements of D is reversed.
+*           Thus if MODE is positive, D has entries ranging from
+*              1 to 1/COND, if negative, from 1/COND to 1,
+*           If SYM='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.
+*
+*  RANK   - INTEGER
+*           The rank of matrix to be generated for modes 1,2,3 only.
+*           D( RANK+1:N ) = 0.
+*           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 ZLATMT 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 DLATM7
+*            2  => Cannot scale to DMAX (max. sing. value is 0)
+*            3  => Error return from ZLAGGE, ZLAGHE or ZLAGSY
+*
+*  =====================================================================
+*
+*     .. 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 ..
+      COMPLEX*16         C, CT, DUMMY, EXTRA, S, ST, ZTEMP
+      DOUBLE PRECISION   ALPHA, ANGLE, REALC, TEMP
+      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
+      LOGICAL            CSYM, GIVENS, ILEXTR, ILTEMP, TOPDWN
+*     ..
+*     .. External Functions ..
+      COMPLEX*16         ZLARND
+      DOUBLE PRECISION   DLARND
+      LOGICAL            LSAME
+      EXTERNAL           ZLARND, DLARND, LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLATM7, 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
+         CSYM = .FALSE.
+      ELSE IF( LSAME( SYM, 'P' ) ) THEN
+         ISYM = 2
+         IRSIGN = 0
+         CSYM = .FALSE.
+      ELSE IF( LSAME( SYM, 'S' ) ) THEN
+         ISYM = 2
+         IRSIGN = 0
+         CSYM = .TRUE.
+      ELSE IF( LSAME( SYM, 'H' ) ) THEN
+         ISYM = 2
+         IRSIGN = 1
+         CSYM = .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 )
+*
+      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( 'ZLATMT', -INFO )
+         RETURN
+      END IF
+*
+*     Initialize random number generator
+*
+      DO 100 I = 1, 4
+         ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
+  100 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 DLATM7( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, RANK,
+     $             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( RANK ) ) ) 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 110 I = 2, RANK
+            TEMP = MAX( TEMP, ABS( D( I ) ) )
+  110    CONTINUE
+*
+         IF( TEMP.GT.ZERO ) THEN
+            ALPHA = DMAX / TEMP
+         ELSE
+            INFO = 2
+            RETURN
+         END IF
+*
+         CALL DSCAL( RANK, 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 120 J = 1, MNMIN
+            A( ( 1-ISKEW )*J+IOFFST, J ) = DCMPLX( D( J ) )
+  120    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 130 J = 1, MNMIN
+               A( ( 1-ISKEW )*J+IOFFST, J ) = DCMPLX( D( J ) )
+  130       CONTINUE
+*
+            IF( TOPDWN ) THEN
+               JKL = 0
+               DO 160 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 150 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 140 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 = DLARND( 5, ISEED )
+                           C = DCONJG( REALC*DUMMY )
+                           S = DCONJG( -S*DUMMY )
+                        END IF
+                        IROW = MAX( 1, JCH-JKU )
+                        IL = IR + 2 - IROW
+                        ZTEMP = CZERO
+                        ILTEMP = JCH.GT.JKU
+                        CALL ZLAROT( .FALSE., ILTEMP, .TRUE., IL, C, S,
+     $                               A( IROW-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, ZTEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL ZLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
+     $                                  IC+1 ), ZTEMP, 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,
+     $                                  ZTEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  140                CONTINUE
+  150             CONTINUE
+  160          CONTINUE
+*
+               JKU = UUB
+               DO 190 JKL = 1, LLB
+*
+*                 Transform from bandwidth JKL-1, JKU to JKL, JKU
+*
+                  DO 180 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 170 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
+                        ZTEMP = CZERO
+                        ILTEMP = JCH.GT.JKL
+                        CALL ZLAROT( .TRUE., ILTEMP, .TRUE., IL, C, S,
+     $                               A( IR-ISKEW*ICOL+IOFFST, ICOL ),
+     $                               ILDA, ZTEMP, EXTRA )
+                        IF( ILTEMP ) THEN
+                           CALL ZLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
+     $                                  ICOL+1 ), ZTEMP, 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,
+     $                                  ZTEMP )
+                           IC = ICOL
+                           IR = IROW
+                        END IF
+  170                CONTINUE
+  180             CONTINUE
+  190          CONTINUE
+*
+            ELSE
+*
+*              Bottom-Up -- Start at the bottom right.
+*
+               JKL = 0
+               DO 220 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 210 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 200 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
+                        ZTEMP = CZERO
+                        CALL ZLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
+     $                               C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
+     $                               ILDA, EXTRA, ZTEMP )
+                        IF( ILTEMP ) THEN
+                           CALL ZLARTG( A( JCH-ISKEW*ICOL+IOFFST,
+     $                                  ICOL ), ZTEMP, 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, ZTEMP, EXTRA )
+                           IC = ICOL
+                        END IF
+  200                CONTINUE
+  210             CONTINUE
+  220          CONTINUE
+*
+               JKU = UUB
+               DO 250 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 240 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 230 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
+                        ZTEMP = CZERO
+                        CALL ZLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
+     $                               C, S, A( IR-ISKEW*JCH+IOFFST,
+     $                               JCH ), ILDA, EXTRA, ZTEMP )
+                        IF( ILTEMP ) THEN
+                           CALL ZLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
+     $                                  ZTEMP, 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, ZTEMP, EXTRA )
+                           IR = IROW
+                        END IF
+  230                CONTINUE
+  240             CONTINUE
+  250          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 260 J = 1, MNMIN
+                  A( ( 1-ISKEW )*J+IOFFG, J ) = DCMPLX( D( J ) )
+  260          CONTINUE
+*
+               DO 290 K = 1, UUB
+                  DO 280 JC = 1, N - 1
+                     IROW = MAX( 1, JC-K )
+                     IL = MIN( JC+1, K+2 )
+                     EXTRA = CZERO
+                     ZTEMP = 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( CSYM ) THEN
+                        CT = C
+                        ST = S
+                     ELSE
+                        ZTEMP = DCONJG( ZTEMP )
+                        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, ZTEMP )
+                     CALL ZLAROT( .TRUE., .TRUE., .FALSE.,
+     $                            MIN( K, N-JC )+1, CT, ST,
+     $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
+     $                            ZTEMP, DUMMY )
+*
+*                    Chase EXTRA back up the matrix
+*
+                     ICOL = JC
+                     DO 270 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 )
+                        ZTEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
+                        IF( CSYM ) THEN
+                           CT = C
+                           ST = S
+                        ELSE
+                           ZTEMP = DCONJG( ZTEMP )
+                           CT = DCONJG( C )
+                           ST = DCONJG( S )
+                        END IF
+                        CALL ZLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
+     $                               A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
+     $                               ILDA, ZTEMP, 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, ZTEMP )
+                        ICOL = JCH
+  270                CONTINUE
+  280             CONTINUE
+  290          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 320 JC = 1, N
+                     IROW = IOFFST - ISKEW*JC
+                     IF( CSYM ) THEN
+                        DO 300 JR = JC, MIN( N, JC+UUB )
+                           A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  300                   CONTINUE
+                     ELSE
+                        DO 310 JR = JC, MIN( N, JC+UUB )
+                           A( JR+IROW, JC ) = DCONJG( A( JC-ISKEW*JR+
+     $                                        IOFFG, JR ) )
+  310                   CONTINUE
+                     END IF
+  320             CONTINUE
+                  IF( IPACK.EQ.5 ) THEN
+                     DO 340 JC = N - UUB + 1, N
+                        DO 330 JR = N + 2 - JC, UUB + 1
+                           A( JR, JC ) = CZERO
+  330                   CONTINUE
+  340                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 350 J = 1, MNMIN
+                  A( ( 1-ISKEW )*J+IOFFG, J ) = DCMPLX( D( J ) )
+  350          CONTINUE
+*
+               DO 380 K = 1, UUB
+                  DO 370 JC = N - 1, 1, -1
+                     IL = MIN( N+1-JC, K+2 )
+                     EXTRA = CZERO
+                     ZTEMP = 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( CSYM ) THEN
+                        CT = C
+                        ST = S
+                     ELSE
+                        ZTEMP = DCONJG( ZTEMP )
+                        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,
+     $                            ZTEMP, 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, ZTEMP )
+*
+*                    Chase EXTRA back down the matrix
+*
+                     ICOL = JC
+                     DO 360 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
+                        ZTEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
+                        IF( CSYM ) THEN
+                           CT = C
+                           ST = S
+                        ELSE
+                           ZTEMP = DCONJG( ZTEMP )
+                           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, ZTEMP )
+                        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, ZTEMP, EXTRA )
+                        ICOL = JCH
+  360                CONTINUE
+  370             CONTINUE
+  380          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 410 JC = N, 1, -1
+                     IROW = IOFFST - ISKEW*JC
+                     IF( CSYM ) THEN
+                        DO 390 JR = JC, MAX( 1, JC-UUB ), -1
+                           A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
+  390                   CONTINUE
+                     ELSE
+                        DO 400 JR = JC, MAX( 1, JC-UUB ), -1
+                           A( JR+IROW, JC ) = DCONJG( A( JC-ISKEW*JR+
+     $                                        IOFFG, JR ) )
+  400                   CONTINUE
+                     END IF
+  410             CONTINUE
+                  IF( IPACK.EQ.6 ) THEN
+                     DO 430 JC = 1, UUB
+                        DO 420 JR = 1, UUB + 1 - JC
+                           A( JR, JC ) = CZERO
+  420                   CONTINUE
+  430                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.CSYM ) THEN
+               DO 440 JC = 1, N
+                  IROW = IOFFST + ( 1-ISKEW )*JC
+                  A( IROW, JC ) = DCMPLX( DBLE( A( IROW, JC ) ) )
+  440          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( CSYM ) 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 460 J = 1, M
+               DO 450 I = J + 1, M
+                  A( I, J ) = CZERO
+  450          CONTINUE
+  460       CONTINUE
+*
+         ELSE IF( IPACK.EQ.2 ) THEN
+*
+*           'L' -- Lower triangular, not packed
+*
+            DO 480 J = 2, M
+               DO 470 I = 1, J - 1
+                  A( I, J ) = CZERO
+  470          CONTINUE
+  480       CONTINUE
+*
+         ELSE IF( IPACK.EQ.3 ) THEN
+*
+*           'C' -- Upper triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 500 J = 1, M
+               DO 490 I = 1, J
+                  IROW = IROW + 1
+                  IF( IROW.GT.LDA ) THEN
+                     IROW = 1
+                     ICOL = ICOL + 1
+                  END IF
+                  A( IROW, ICOL ) = A( I, J )
+  490          CONTINUE
+  500       CONTINUE
+*
+         ELSE IF( IPACK.EQ.4 ) THEN
+*
+*           'R' -- Lower triangle packed Columnwise.
+*
+            ICOL = 1
+            IROW = 0
+            DO 520 J = 1, M
+               DO 510 I = J, M
+                  IROW = IROW + 1
+                  IF( IROW.GT.LDA ) THEN
+                     IROW = 1
+                     ICOL = ICOL + 1
+                  END IF
+                  A( IROW, ICOL ) = A( I, J )
+  510          CONTINUE
+  520       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 540 J = 1, UUB
+               DO 530 I = MIN( J+LLB, M ), 1, -1
+                  A( I-J+UUB+1, J ) = A( I, J )
+  530          CONTINUE
+  540       CONTINUE
+*
+            DO 560 J = UUB + 2, N
+               DO 550 I = J - UUB, MIN( J+LLB, M )
+                  A( I-J+UUB+1, J ) = A( I, J )
+  550          CONTINUE
+  560       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 580 JC = ICOL, M
+               DO 570 JR = IROW + 1, LDA
+                  A( JR, JC ) = CZERO
+  570          CONTINUE
+               IROW = 0
+  580       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 610 JC = 1, N
+               DO 590 JR = 1, UUB + 1 - JC
+                  A( JR, JC ) = CZERO
+  590          CONTINUE
+               DO 600 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
+                  A( JR, JC ) = CZERO
+  600          CONTINUE
+  610       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of ZLATMT
+*
+      END