LAPACK: TESTING/MATGEN/dlatms.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
00002      $                   KL, KU, PACK, A, LDA, WORK, INFO )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     June 2010
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          DIST, PACK, SYM
00010       INTEGER            INFO, KL, KU, LDA, M, MODE, N
00011       DOUBLE PRECISION   COND, DMAX
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            ISEED( 4 )
00015       DOUBLE PRECISION   A( LDA, * ), D( * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *     DLATMS generates random matrices with specified singular values
00022 *     (or symmetric/hermitian with specified eigenvalues)
00023 *     for testing LAPACK programs.
00024 *
00025 *     DLATMS operates by applying the following sequence of
00026 *     operations:
00027 *
00028 *       Set the diagonal to D, where D may be input or
00029 *          computed according to MODE, COND, DMAX, and SYM
00030 *          as described below.
00031 *
00032 *       Generate a matrix with the appropriate band structure, by one
00033 *          of two methods:
00034 *
00035 *       Method A:
00036 *           Generate a dense M x N matrix by multiplying D on the left
00037 *               and the right by random unitary matrices, then:
00038 *
00039 *           Reduce the bandwidth according to KL and KU, using
00040 *           Householder transformations.
00041 *
00042 *       Method B:
00043 *           Convert the bandwidth-0 (i.e., diagonal) matrix to a
00044 *               bandwidth-1 matrix using Givens rotations, "chasing"
00045 *               out-of-band elements back, much as in QR; then
00046 *               convert the bandwidth-1 to a bandwidth-2 matrix, etc.
00047 *               Note that for reasonably small bandwidths (relative to
00048 *               M and N) this requires less storage, as a dense matrix
00049 *               is not generated.  Also, for symmetric matrices, only
00050 *               one triangle is generated.
00051 *
00052 *       Method A is chosen if the bandwidth is a large fraction of the
00053 *           order of the matrix, and LDA is at least M (so a dense
00054 *           matrix can be stored.)  Method B is chosen if the bandwidth
00055 *           is small (< 1/2 N for symmetric, < .3 N+M for
00056 *           non-symmetric), or LDA is less than M and not less than the
00057 *           bandwidth.
00058 *
00059 *       Pack the matrix if desired. Options specified by PACK are:
00060 *          no packing
00061 *          zero out upper half (if symmetric)
00062 *          zero out lower half (if symmetric)
00063 *          store the upper half columnwise (if symmetric or upper
00064 *                triangular)
00065 *          store the lower half columnwise (if symmetric or lower
00066 *                triangular)
00067 *          store the lower triangle in banded format (if symmetric
00068 *                or lower triangular)
00069 *          store the upper triangle in banded format (if symmetric
00070 *                or upper triangular)
00071 *          store the entire matrix in banded format
00072 *       If Method B is chosen, and band format is specified, then the
00073 *          matrix will be generated in the band format, so no repacking
00074 *          will be necessary.
00075 *
00076 *  Arguments
00077 *  =========
00078 *
00079 *  M        (input) INTEGER
00080 *           The number of rows of A. Not modified.
00081 *
00082 *  N        (input) INTEGER
00083 *           The number of columns of A. Not modified.
00084 *
00085 *  DIST     (input) CHARACTER*1
00086 *           On entry, DIST specifies the type of distribution to be used
00087 *           to generate the random eigen-/singular values.
00088 *           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
00089 *           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
00090 *           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
00091 *           Not modified.
00092 *
00093 *  ISEED    (input/output) INTEGER array, dimension ( 4 )
00094 *           On entry ISEED specifies the seed of the random number
00095 *           generator. They should lie between 0 and 4095 inclusive,
00096 *           and ISEED(4) should be odd. The random number generator
00097 *           uses a linear congruential sequence limited to small
00098 *           integers, and so should produce machine independent
00099 *           random numbers. The values of ISEED are changed on
00100 *           exit, and can be used in the next call to DLATMS
00101 *           to continue the same random number sequence.
00102 *           Changed on exit.
00103 *
00104 *  SYM      (input) CHARACTER*1
00105 *           If SYM='S' or 'H', the generated matrix is symmetric, with
00106 *             eigenvalues specified by D, COND, MODE, and DMAX; they
00107 *             may be positive, negative, or zero.
00108 *           If SYM='P', the generated matrix is symmetric, with
00109 *             eigenvalues (= singular values) specified by D, COND,
00110 *             MODE, and DMAX; they will not be negative.
00111 *           If SYM='N', the generated matrix is nonsymmetric, with
00112 *             singular values specified by D, COND, MODE, and DMAX;
00113 *             they will not be negative.
00114 *           Not modified.
00115 *
00116 *  D        (input/output) DOUBLE PRECISION array, dimension ( MIN( M , N ) )
00117 *           This array is used to specify the singular values or
00118 *           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
00119 *           assumed to contain the singular/eigenvalues, otherwise
00120 *           they will be computed according to MODE, COND, and DMAX,
00121 *           and placed in D.
00122 *           Modified if MODE is nonzero.
00123 *
00124 *  MODE     (input) INTEGER
00125 *           On entry this describes how the singular/eigenvalues are to
00126 *           be specified:
00127 *           MODE = 0 means use D as input
00128 *           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
00129 *           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
00130 *           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
00131 *           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
00132 *           MODE = 5 sets D to random numbers in the range
00133 *                    ( 1/COND , 1 ) such that their logarithms
00134 *                    are uniformly distributed.
00135 *           MODE = 6 set D to random numbers from same distribution
00136 *                    as the rest of the matrix.
00137 *           MODE < 0 has the same meaning as ABS(MODE), except that
00138 *              the order of the elements of D is reversed.
00139 *           Thus if MODE is positive, D has entries ranging from
00140 *              1 to 1/COND, if negative, from 1/COND to 1,
00141 *           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
00142 *              the elements of D will also be multiplied by a random
00143 *              sign (i.e., +1 or -1.)
00144 *           Not modified.
00145 *
00146 *  COND     (input) DOUBLE PRECISION
00147 *           On entry, this is used as described under MODE above.
00148 *           If used, it must be >= 1. Not modified.
00149 *
00150 *  DMAX     (input) DOUBLE PRECISION
00151 *           If MODE is neither -6, 0 nor 6, the contents of D, as
00152 *           computed according to MODE and COND, will be scaled by
00153 *           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
00154 *           singular value (which is to say the norm) will be abs(DMAX).
00155 *           Note that DMAX need not be positive: if DMAX is negative
00156 *           (or zero), D will be scaled by a negative number (or zero).
00157 *           Not modified.
00158 *
00159 *  KL       (input) INTEGER
00160 *           This specifies the lower bandwidth of the  matrix. For
00161 *           example, KL=0 implies upper triangular, KL=1 implies upper
00162 *           Hessenberg, and KL being at least M-1 means that the matrix
00163 *           has full lower bandwidth.  KL must equal KU if the matrix
00164 *           is symmetric.
00165 *           Not modified.
00166 *
00167 *  KU       (input) INTEGER
00168 *           This specifies the upper bandwidth of the  matrix. For
00169 *           example, KU=0 implies lower triangular, KU=1 implies lower
00170 *           Hessenberg, and KU being at least N-1 means that the matrix
00171 *           has full upper bandwidth.  KL must equal KU if the matrix
00172 *           is symmetric.
00173 *           Not modified.
00174 *
00175 *  PACK     (input) CHARACTER*1
00176 *           This specifies packing of matrix as follows:
00177 *           'N' => no packing
00178 *           'U' => zero out all subdiagonal entries (if symmetric)
00179 *           'L' => zero out all superdiagonal entries (if symmetric)
00180 *           'C' => store the upper triangle columnwise
00181 *                  (only if the matrix is symmetric or upper triangular)
00182 *           'R' => store the lower triangle columnwise
00183 *                  (only if the matrix is symmetric or lower triangular)
00184 *           'B' => store the lower triangle in band storage scheme
00185 *                  (only if matrix symmetric or lower triangular)
00186 *           'Q' => store the upper triangle in band storage scheme
00187 *                  (only if matrix symmetric or upper triangular)
00188 *           'Z' => store the entire matrix in band storage scheme
00189 *                      (pivoting can be provided for by using this
00190 *                      option to store A in the trailing rows of
00191 *                      the allocated storage)
00192 *
00193 *           Using these options, the various LAPACK packed and banded
00194 *           storage schemes can be obtained:
00195 *           GB               - use 'Z'
00196 *           PB, SB or TB     - use 'B' or 'Q'
00197 *           PP, SP or TP     - use 'C' or 'R'
00198 *
00199 *           If two calls to DLATMS differ only in the PACK parameter,
00200 *           they will generate mathematically equivalent matrices.
00201 *           Not modified.
00202 *
00203 *  A        (input/output) DOUBLE PRECISION array, dimension ( LDA, N )
00204 *           On exit A is the desired test matrix.  A is first generated
00205 *           in full (unpacked) form, and then packed, if so specified
00206 *           by PACK.  Thus, the first M elements of the first N
00207 *           columns will always be modified.  If PACK specifies a
00208 *           packed or banded storage scheme, all LDA elements of the
00209 *           first N columns will be modified; the elements of the
00210 *           array which do not correspond to elements of the generated
00211 *           matrix are set to zero.
00212 *           Modified.
00213 *
00214 *  LDA      (input) INTEGER
00215 *           LDA specifies the first dimension of A as declared in the
00216 *           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
00217 *           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
00218 *           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
00219 *           If PACK='Z', LDA must be large enough to hold the packed
00220 *           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
00221 *           Not modified.
00222 *
00223 *  WORK     (workspace) DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) )
00224 *           Workspace.
00225 *           Modified.
00226 *
00227 *  INFO     (output) INTEGER
00228 *           Error code.  On exit, INFO will be set to one of the
00229 *           following values:
00230 *             0 => normal return
00231 *            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
00232 *            -2 => N negative
00233 *            -3 => DIST illegal string
00234 *            -5 => SYM illegal string
00235 *            -7 => MODE not in range -6 to 6
00236 *            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
00237 *           -10 => KL negative
00238 *           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
00239 *           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
00240 *                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
00241 *                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
00242 *                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
00243 *                  N.
00244 *           -14 => LDA is less than M, or PACK='Z' and LDA is less than
00245 *                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
00246 *            1  => Error return from DLATM1
00247 *            2  => Cannot scale to DMAX (max. sing. value is 0)
00248 *            3  => Error return from DLAGGE or SLAGSY
00249 *
00250 *  =====================================================================
00251 *
00252 *     .. Parameters ..
00253       DOUBLE PRECISION   ZERO
00254       PARAMETER          ( ZERO = 0.0D0 )
00255       DOUBLE PRECISION   ONE
00256       PARAMETER          ( ONE = 1.0D0 )
00257       DOUBLE PRECISION   TWOPI
00258       PARAMETER          ( TWOPI = 6.2831853071795864769252867663D+0 )
00259 *     ..
00260 *     .. Local Scalars ..
00261       LOGICAL            GIVENS, ILEXTR, ILTEMP, TOPDWN
00262       INTEGER            I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
00263      $                   IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2,
00264      $                   IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH,
00265      $                   JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC,
00266      $                   UUB
00267       DOUBLE PRECISION   ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
00268 *     ..
00269 *     .. External Functions ..
00270       LOGICAL            LSAME
00271       DOUBLE PRECISION   DLARND
00272       EXTERNAL           LSAME, DLARND
00273 *     ..
00274 *     .. External Subroutines ..
00275       EXTERNAL           DCOPY, DLAGGE, DLAGSY, DLAROT, DLARTG, DLASET,
00276      $                   DLATM1, DSCAL, XERBLA
00277 *     ..
00278 *     .. Intrinsic Functions ..
00279       INTRINSIC          ABS, COS, DBLE, MAX, MIN, MOD, SIN
00280 *     ..
00281 *     .. Executable Statements ..
00282 *
00283 *     1)      Decode and Test the input parameters.
00284 *             Initialize flags & seed.
00285 *
00286       INFO = 0
00287 *
00288 *     Quick return if possible
00289 *
00290       IF( M.EQ.0 .OR. N.EQ.0 )
00291      $   RETURN
00292 *
00293 *     Decode DIST
00294 *
00295       IF( LSAME( DIST, 'U' ) ) THEN
00296          IDIST = 1
00297       ELSE IF( LSAME( DIST, 'S' ) ) THEN
00298          IDIST = 2
00299       ELSE IF( LSAME( DIST, 'N' ) ) THEN
00300          IDIST = 3
00301       ELSE
00302          IDIST = -1
00303       END IF
00304 *
00305 *     Decode SYM
00306 *
00307       IF( LSAME( SYM, 'N' ) ) THEN
00308          ISYM = 1
00309          IRSIGN = 0
00310       ELSE IF( LSAME( SYM, 'P' ) ) THEN
00311          ISYM = 2
00312          IRSIGN = 0
00313       ELSE IF( LSAME( SYM, 'S' ) ) THEN
00314          ISYM = 2
00315          IRSIGN = 1
00316       ELSE IF( LSAME( SYM, 'H' ) ) THEN
00317          ISYM = 2
00318          IRSIGN = 1
00319       ELSE
00320          ISYM = -1
00321       END IF
00322 *
00323 *     Decode PACK
00324 *
00325       ISYMPK = 0
00326       IF( LSAME( PACK, 'N' ) ) THEN
00327          IPACK = 0
00328       ELSE IF( LSAME( PACK, 'U' ) ) THEN
00329          IPACK = 1
00330          ISYMPK = 1
00331       ELSE IF( LSAME( PACK, 'L' ) ) THEN
00332          IPACK = 2
00333          ISYMPK = 1
00334       ELSE IF( LSAME( PACK, 'C' ) ) THEN
00335          IPACK = 3
00336          ISYMPK = 2
00337       ELSE IF( LSAME( PACK, 'R' ) ) THEN
00338          IPACK = 4
00339          ISYMPK = 3
00340       ELSE IF( LSAME( PACK, 'B' ) ) THEN
00341          IPACK = 5
00342          ISYMPK = 3
00343       ELSE IF( LSAME( PACK, 'Q' ) ) THEN
00344          IPACK = 6
00345          ISYMPK = 2
00346       ELSE IF( LSAME( PACK, 'Z' ) ) THEN
00347          IPACK = 7
00348       ELSE
00349          IPACK = -1
00350       END IF
00351 *
00352 *     Set certain internal parameters
00353 *
00354       MNMIN = MIN( M, N )
00355       LLB = MIN( KL, M-1 )
00356       UUB = MIN( KU, N-1 )
00357       MR = MIN( M, N+LLB )
00358       NC = MIN( N, M+UUB )
00359 *
00360       IF( IPACK.EQ.5 .OR. IPACK.EQ.6 ) THEN
00361          MINLDA = UUB + 1
00362       ELSE IF( IPACK.EQ.7 ) THEN
00363          MINLDA = LLB + UUB + 1
00364       ELSE
00365          MINLDA = M
00366       END IF
00367 *
00368 *     Use Givens rotation method if bandwidth small enough,
00369 *     or if LDA is too small to store the matrix unpacked.
00370 *
00371       GIVENS = .FALSE.
00372       IF( ISYM.EQ.1 ) THEN
00373          IF( DBLE( LLB+UUB ).LT.0.3D0*DBLE( MAX( 1, MR+NC ) ) )
00374      $      GIVENS = .TRUE.
00375       ELSE
00376          IF( 2*LLB.LT.M )
00377      $      GIVENS = .TRUE.
00378       END IF
00379       IF( LDA.LT.M .AND. LDA.GE.MINLDA )
00380      $   GIVENS = .TRUE.
00381 *
00382 *     Set INFO if an error
00383 *
00384       IF( M.LT.0 ) THEN
00385          INFO = -1
00386       ELSE IF( M.NE.N .AND. ISYM.NE.1 ) THEN
00387          INFO = -1
00388       ELSE IF( N.LT.0 ) THEN
00389          INFO = -2
00390       ELSE IF( IDIST.EQ.-1 ) THEN
00391          INFO = -3
00392       ELSE IF( ISYM.EQ.-1 ) THEN
00393          INFO = -5
00394       ELSE IF( ABS( MODE ).GT.6 ) THEN
00395          INFO = -7
00396       ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
00397      $          THEN
00398          INFO = -8
00399       ELSE IF( KL.LT.0 ) THEN
00400          INFO = -10
00401       ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN
00402          INFO = -11
00403       ELSE IF( IPACK.EQ.-1 .OR. ( ISYMPK.EQ.1 .AND. ISYM.EQ.1 ) .OR.
00404      $         ( ISYMPK.EQ.2 .AND. ISYM.EQ.1 .AND. KL.GT.0 ) .OR.
00405      $         ( ISYMPK.EQ.3 .AND. ISYM.EQ.1 .AND. KU.GT.0 ) .OR.
00406      $         ( ISYMPK.NE.0 .AND. M.NE.N ) ) THEN
00407          INFO = -12
00408       ELSE IF( LDA.LT.MAX( 1, MINLDA ) ) THEN
00409          INFO = -14
00410       END IF
00411 *
00412       IF( INFO.NE.0 ) THEN
00413          CALL XERBLA( 'DLATMS', -INFO )
00414          RETURN
00415       END IF
00416 *
00417 *     Initialize random number generator
00418 *
00419       DO 10 I = 1, 4
00420          ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
00421    10 CONTINUE
00422 *
00423       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
00424      $   ISEED( 4 ) = ISEED( 4 ) + 1
00425 *
00426 *     2)      Set up D  if indicated.
00427 *
00428 *             Compute D according to COND and MODE
00429 *
00430       CALL DLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, IINFO )
00431       IF( IINFO.NE.0 ) THEN
00432          INFO = 1
00433          RETURN
00434       END IF
00435 *
00436 *     Choose Top-Down if D is (apparently) increasing,
00437 *     Bottom-Up if D is (apparently) decreasing.
00438 *
00439       IF( ABS( D( 1 ) ).LE.ABS( D( MNMIN ) ) ) THEN
00440          TOPDWN = .TRUE.
00441       ELSE
00442          TOPDWN = .FALSE.
00443       END IF
00444 *
00445       IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
00446 *
00447 *        Scale by DMAX
00448 *
00449          TEMP = ABS( D( 1 ) )
00450          DO 20 I = 2, MNMIN
00451             TEMP = MAX( TEMP, ABS( D( I ) ) )
00452    20    CONTINUE
00453 *
00454          IF( TEMP.GT.ZERO ) THEN
00455             ALPHA = DMAX / TEMP
00456          ELSE
00457             INFO = 2
00458             RETURN
00459          END IF
00460 *
00461          CALL DSCAL( MNMIN, ALPHA, D, 1 )
00462 *
00463       END IF
00464 *
00465 *     3)      Generate Banded Matrix using Givens rotations.
00466 *             Also the special case of UUB=LLB=0
00467 *
00468 *               Compute Addressing constants to cover all
00469 *               storage formats.  Whether GE, SY, GB, or SB,
00470 *               upper or lower triangle or both,
00471 *               the (i,j)-th element is in
00472 *               A( i - ISKEW*j + IOFFST, j )
00473 *
00474       IF( IPACK.GT.4 ) THEN
00475          ILDA = LDA - 1
00476          ISKEW = 1
00477          IF( IPACK.GT.5 ) THEN
00478             IOFFST = UUB + 1
00479          ELSE
00480             IOFFST = 1
00481          END IF
00482       ELSE
00483          ILDA = LDA
00484          ISKEW = 0
00485          IOFFST = 0
00486       END IF
00487 *
00488 *     IPACKG is the format that the matrix is generated in. If this is
00489 *     different from IPACK, then the matrix must be repacked at the
00490 *     end.  It also signals how to compute the norm, for scaling.
00491 *
00492       IPACKG = 0
00493       CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
00494 *
00495 *     Diagonal Matrix -- We are done, unless it
00496 *     is to be stored SP/PP/TP (PACK='R' or 'C')
00497 *
00498       IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN
00499          CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
00500          IF( IPACK.LE.2 .OR. IPACK.GE.5 )
00501      $      IPACKG = IPACK
00502 *
00503       ELSE IF( GIVENS ) THEN
00504 *
00505 *        Check whether to use Givens rotations,
00506 *        Householder transformations, or nothing.
00507 *
00508          IF( ISYM.EQ.1 ) THEN
00509 *
00510 *           Non-symmetric -- A = U D V
00511 *
00512             IF( IPACK.GT.4 ) THEN
00513                IPACKG = IPACK
00514             ELSE
00515                IPACKG = 0
00516             END IF
00517 *
00518             CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
00519 *
00520             IF( TOPDWN ) THEN
00521                JKL = 0
00522                DO 50 JKU = 1, UUB
00523 *
00524 *                 Transform from bandwidth JKL, JKU-1 to JKL, JKU
00525 *
00526 *                 Last row actually rotated is M
00527 *                 Last column actually rotated is MIN( M+JKU, N )
00528 *
00529                   DO 40 JR = 1, MIN( M+JKU, N ) + JKL - 1
00530                      EXTRA = ZERO
00531                      ANGLE = TWOPI*DLARND( 1, ISEED )
00532                      C = COS( ANGLE )
00533                      S = SIN( ANGLE )
00534                      ICOL = MAX( 1, JR-JKL )
00535                      IF( JR.LT.M ) THEN
00536                         IL = MIN( N, JR+JKU ) + 1 - ICOL
00537                         CALL DLAROT( .TRUE., JR.GT.JKL, .FALSE., IL, C,
00538      $                               S, A( JR-ISKEW*ICOL+IOFFST, ICOL ),
00539      $                               ILDA, EXTRA, DUMMY )
00540                      END IF
00541 *
00542 *                    Chase "EXTRA" back up
00543 *
00544                      IR = JR
00545                      IC = ICOL
00546                      DO 30 JCH = JR - JKL, 1, -JKL - JKU
00547                         IF( IR.LT.M ) THEN
00548                            CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
00549      $                                  IC+1 ), EXTRA, C, S, DUMMY )
00550                         END IF
00551                         IROW = MAX( 1, JCH-JKU )
00552                         IL = IR + 2 - IROW
00553                         TEMP = ZERO
00554                         ILTEMP = JCH.GT.JKU
00555                         CALL DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S,
00556      $                               A( IROW-ISKEW*IC+IOFFST, IC ),
00557      $                               ILDA, TEMP, EXTRA )
00558                         IF( ILTEMP ) THEN
00559                            CALL DLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
00560      $                                  IC+1 ), TEMP, C, S, DUMMY )
00561                            ICOL = MAX( 1, JCH-JKU-JKL )
00562                            IL = IC + 2 - ICOL
00563                            EXTRA = ZERO
00564                            CALL DLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
00565      $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
00566      $                                  IOFFST, ICOL ), ILDA, EXTRA,
00567      $                                  TEMP )
00568                            IC = ICOL
00569                            IR = IROW
00570                         END IF
00571    30                CONTINUE
00572    40             CONTINUE
00573    50          CONTINUE
00574 *
00575                JKU = UUB
00576                DO 80 JKL = 1, LLB
00577 *
00578 *                 Transform from bandwidth JKL-1, JKU to JKL, JKU
00579 *
00580                   DO 70 JC = 1, MIN( N+JKL, M ) + JKU - 1
00581                      EXTRA = ZERO
00582                      ANGLE = TWOPI*DLARND( 1, ISEED )
00583                      C = COS( ANGLE )
00584                      S = SIN( ANGLE )
00585                      IROW = MAX( 1, JC-JKU )
00586                      IF( JC.LT.N ) THEN
00587                         IL = MIN( M, JC+JKL ) + 1 - IROW
00588                         CALL DLAROT( .FALSE., JC.GT.JKU, .FALSE., IL, C,
00589      $                               S, A( IROW-ISKEW*JC+IOFFST, JC ),
00590      $                               ILDA, EXTRA, DUMMY )
00591                      END IF
00592 *
00593 *                    Chase "EXTRA" back up
00594 *
00595                      IC = JC
00596                      IR = IROW
00597                      DO 60 JCH = JC - JKU, 1, -JKL - JKU
00598                         IF( IC.LT.N ) THEN
00599                            CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
00600      $                                  IC+1 ), EXTRA, C, S, DUMMY )
00601                         END IF
00602                         ICOL = MAX( 1, JCH-JKL )
00603                         IL = IC + 2 - ICOL
00604                         TEMP = ZERO
00605                         ILTEMP = JCH.GT.JKL
00606                         CALL DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S,
00607      $                               A( IR-ISKEW*ICOL+IOFFST, ICOL ),
00608      $                               ILDA, TEMP, EXTRA )
00609                         IF( ILTEMP ) THEN
00610                            CALL DLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
00611      $                                  ICOL+1 ), TEMP, C, S, DUMMY )
00612                            IROW = MAX( 1, JCH-JKL-JKU )
00613                            IL = IR + 2 - IROW
00614                            EXTRA = ZERO
00615                            CALL DLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
00616      $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
00617      $                                  IOFFST, ICOL ), ILDA, EXTRA,
00618      $                                  TEMP )
00619                            IC = ICOL
00620                            IR = IROW
00621                         END IF
00622    60                CONTINUE
00623    70             CONTINUE
00624    80          CONTINUE
00625 *
00626             ELSE
00627 *
00628 *              Bottom-Up -- Start at the bottom right.
00629 *
00630                JKL = 0
00631                DO 110 JKU = 1, UUB
00632 *
00633 *                 Transform from bandwidth JKL, JKU-1 to JKL, JKU
00634 *
00635 *                 First row actually rotated is M
00636 *                 First column actually rotated is MIN( M+JKU, N )
00637 *
00638                   IENDCH = MIN( M, N+JKL ) - 1
00639                   DO 100 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
00640                      EXTRA = ZERO
00641                      ANGLE = TWOPI*DLARND( 1, ISEED )
00642                      C = COS( ANGLE )
00643                      S = SIN( ANGLE )
00644                      IROW = MAX( 1, JC-JKU+1 )
00645                      IF( JC.GT.0 ) THEN
00646                         IL = MIN( M, JC+JKL+1 ) + 1 - IROW
00647                         CALL DLAROT( .FALSE., .FALSE., JC+JKL.LT.M, IL,
00648      $                               C, S, A( IROW-ISKEW*JC+IOFFST,
00649      $                               JC ), ILDA, DUMMY, EXTRA )
00650                      END IF
00651 *
00652 *                    Chase "EXTRA" back down
00653 *
00654                      IC = JC
00655                      DO 90 JCH = JC + JKL, IENDCH, JKL + JKU
00656                         ILEXTR = IC.GT.0
00657                         IF( ILEXTR ) THEN
00658                            CALL DLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
00659      $                                  EXTRA, C, S, DUMMY )
00660                         END IF
00661                         IC = MAX( 1, IC )
00662                         ICOL = MIN( N-1, JCH+JKU )
00663                         ILTEMP = JCH + JKU.LT.N
00664                         TEMP = ZERO
00665                         CALL DLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
00666      $                               C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
00667      $                               ILDA, EXTRA, TEMP )
00668                         IF( ILTEMP ) THEN
00669                            CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFST,
00670      $                                  ICOL ), TEMP, C, S, DUMMY )
00671                            IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
00672                            EXTRA = ZERO
00673                            CALL DLAROT( .FALSE., .TRUE.,
00674      $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
00675      $                                  A( JCH-ISKEW*ICOL+IOFFST,
00676      $                                  ICOL ), ILDA, TEMP, EXTRA )
00677                            IC = ICOL
00678                         END IF
00679    90                CONTINUE
00680   100             CONTINUE
00681   110          CONTINUE
00682 *
00683                JKU = UUB
00684                DO 140 JKL = 1, LLB
00685 *
00686 *                 Transform from bandwidth JKL-1, JKU to JKL, JKU
00687 *
00688 *                 First row actually rotated is MIN( N+JKL, M )
00689 *                 First column actually rotated is N
00690 *
00691                   IENDCH = MIN( N, M+JKU ) - 1
00692                   DO 130 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
00693                      EXTRA = ZERO
00694                      ANGLE = TWOPI*DLARND( 1, ISEED )
00695                      C = COS( ANGLE )
00696                      S = SIN( ANGLE )
00697                      ICOL = MAX( 1, JR-JKL+1 )
00698                      IF( JR.GT.0 ) THEN
00699                         IL = MIN( N, JR+JKU+1 ) + 1 - ICOL
00700                         CALL DLAROT( .TRUE., .FALSE., JR+JKU.LT.N, IL,
00701      $                               C, S, A( JR-ISKEW*ICOL+IOFFST,
00702      $                               ICOL ), ILDA, DUMMY, EXTRA )
00703                      END IF
00704 *
00705 *                    Chase "EXTRA" back down
00706 *
00707                      IR = JR
00708                      DO 120 JCH = JR + JKU, IENDCH, JKL + JKU
00709                         ILEXTR = IR.GT.0
00710                         IF( ILEXTR ) THEN
00711                            CALL DLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ),
00712      $                                  EXTRA, C, S, DUMMY )
00713                         END IF
00714                         IR = MAX( 1, IR )
00715                         IROW = MIN( M-1, JCH+JKL )
00716                         ILTEMP = JCH + JKL.LT.M
00717                         TEMP = ZERO
00718                         CALL DLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
00719      $                               C, S, A( IR-ISKEW*JCH+IOFFST,
00720      $                               JCH ), ILDA, EXTRA, TEMP )
00721                         IF( ILTEMP ) THEN
00722                            CALL DLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
00723      $                                  TEMP, C, S, DUMMY )
00724                            IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
00725                            EXTRA = ZERO
00726                            CALL DLAROT( .TRUE., .TRUE.,
00727      $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
00728      $                                  A( IROW-ISKEW*JCH+IOFFST, JCH ),
00729      $                                  ILDA, TEMP, EXTRA )
00730                            IR = IROW
00731                         END IF
00732   120                CONTINUE
00733   130             CONTINUE
00734   140          CONTINUE
00735             END IF
00736 *
00737          ELSE
00738 *
00739 *           Symmetric -- A = U D U'
00740 *
00741             IPACKG = IPACK
00742             IOFFG = IOFFST
00743 *
00744             IF( TOPDWN ) THEN
00745 *
00746 *              Top-Down -- Generate Upper triangle only
00747 *
00748                IF( IPACK.GE.5 ) THEN
00749                   IPACKG = 6
00750                   IOFFG = UUB + 1
00751                ELSE
00752                   IPACKG = 1
00753                END IF
00754                CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
00755 *
00756                DO 170 K = 1, UUB
00757                   DO 160 JC = 1, N - 1
00758                      IROW = MAX( 1, JC-K )
00759                      IL = MIN( JC+1, K+2 )
00760                      EXTRA = ZERO
00761                      TEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 )
00762                      ANGLE = TWOPI*DLARND( 1, ISEED )
00763                      C = COS( ANGLE )
00764                      S = SIN( ANGLE )
00765                      CALL DLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
00766      $                            A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
00767      $                            EXTRA, TEMP )
00768                      CALL DLAROT( .TRUE., .TRUE., .FALSE.,
00769      $                            MIN( K, N-JC )+1, C, S,
00770      $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
00771      $                            TEMP, DUMMY )
00772 *
00773 *                    Chase EXTRA back up the matrix
00774 *
00775                      ICOL = JC
00776                      DO 150 JCH = JC - K, 1, -K
00777                         CALL DLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
00778      $                               ICOL+1 ), EXTRA, C, S, DUMMY )
00779                         TEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
00780                         CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, -S,
00781      $                               A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
00782      $                               ILDA, TEMP, EXTRA )
00783                         IROW = MAX( 1, JCH-K )
00784                         IL = MIN( JCH+1, K+2 )
00785                         EXTRA = ZERO
00786                         CALL DLAROT( .FALSE., JCH.GT.K, .TRUE., IL, C,
00787      $                               -S, A( IROW-ISKEW*JCH+IOFFG, JCH ),
00788      $                               ILDA, EXTRA, TEMP )
00789                         ICOL = JCH
00790   150                CONTINUE
00791   160             CONTINUE
00792   170          CONTINUE
00793 *
00794 *              If we need lower triangle, copy from upper. Note that
00795 *              the order of copying is chosen to work for 'q' -> 'b'
00796 *
00797                IF( IPACK.NE.IPACKG .AND. IPACK.NE.3 ) THEN
00798                   DO 190 JC = 1, N
00799                      IROW = IOFFST - ISKEW*JC
00800                      DO 180 JR = JC, MIN( N, JC+UUB )
00801                         A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
00802   180                CONTINUE
00803   190             CONTINUE
00804                   IF( IPACK.EQ.5 ) THEN
00805                      DO 210 JC = N - UUB + 1, N
00806                         DO 200 JR = N + 2 - JC, UUB + 1
00807                            A( JR, JC ) = ZERO
00808   200                   CONTINUE
00809   210                CONTINUE
00810                   END IF
00811                   IF( IPACKG.EQ.6 ) THEN
00812                      IPACKG = IPACK
00813                   ELSE
00814                      IPACKG = 0
00815                   END IF
00816                END IF
00817             ELSE
00818 *
00819 *              Bottom-Up -- Generate Lower triangle only
00820 *
00821                IF( IPACK.GE.5 ) THEN
00822                   IPACKG = 5
00823                   IF( IPACK.EQ.6 )
00824      $               IOFFG = 1
00825                ELSE
00826                   IPACKG = 2
00827                END IF
00828                CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
00829 *
00830                DO 240 K = 1, UUB
00831                   DO 230 JC = N - 1, 1, -1
00832                      IL = MIN( N+1-JC, K+2 )
00833                      EXTRA = ZERO
00834                      TEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
00835                      ANGLE = TWOPI*DLARND( 1, ISEED )
00836                      C = COS( ANGLE )
00837                      S = -SIN( ANGLE )
00838                      CALL DLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S,
00839      $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
00840      $                            TEMP, EXTRA )
00841                      ICOL = MAX( 1, JC-K+1 )
00842                      CALL DLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, C,
00843      $                            S, A( JC-ISKEW*ICOL+IOFFG, ICOL ),
00844      $                            ILDA, DUMMY, TEMP )
00845 *
00846 *                    Chase EXTRA back down the matrix
00847 *
00848                      ICOL = JC
00849                      DO 220 JCH = JC + K, N - 1, K
00850                         CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
00851      $                               EXTRA, C, S, DUMMY )
00852                         TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
00853                         CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
00854      $                               A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
00855      $                               ILDA, EXTRA, TEMP )
00856                         IL = MIN( N+1-JCH, K+2 )
00857                         EXTRA = ZERO
00858                         CALL DLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, C,
00859      $                               S, A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
00860      $                               ILDA, TEMP, EXTRA )
00861                         ICOL = JCH
00862   220                CONTINUE
00863   230             CONTINUE
00864   240          CONTINUE
00865 *
00866 *              If we need upper triangle, copy from lower. Note that
00867 *              the order of copying is chosen to work for 'b' -> 'q'
00868 *
00869                IF( IPACK.NE.IPACKG .AND. IPACK.NE.4 ) THEN
00870                   DO 260 JC = N, 1, -1
00871                      IROW = IOFFST - ISKEW*JC
00872                      DO 250 JR = JC, MAX( 1, JC-UUB ), -1
00873                         A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
00874   250                CONTINUE
00875   260             CONTINUE
00876                   IF( IPACK.EQ.6 ) THEN
00877                      DO 280 JC = 1, UUB
00878                         DO 270 JR = 1, UUB + 1 - JC
00879                            A( JR, JC ) = ZERO
00880   270                   CONTINUE
00881   280                CONTINUE
00882                   END IF
00883                   IF( IPACKG.EQ.5 ) THEN
00884                      IPACKG = IPACK
00885                   ELSE
00886                      IPACKG = 0
00887                   END IF
00888                END IF
00889             END IF
00890          END IF
00891 *
00892       ELSE
00893 *
00894 *        4)      Generate Banded Matrix by first
00895 *                Rotating by random Unitary matrices,
00896 *                then reducing the bandwidth using Householder
00897 *                transformations.
00898 *
00899 *                Note: we should get here only if LDA .ge. N
00900 *
00901          IF( ISYM.EQ.1 ) THEN
00902 *
00903 *           Non-symmetric -- A = U D V
00904 *
00905             CALL DLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
00906      $                   IINFO )
00907          ELSE
00908 *
00909 *           Symmetric -- A = U D U'
00910 *
00911             CALL DLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
00912 *
00913          END IF
00914          IF( IINFO.NE.0 ) THEN
00915             INFO = 3
00916             RETURN
00917          END IF
00918       END IF
00919 *
00920 *     5)      Pack the matrix
00921 *
00922       IF( IPACK.NE.IPACKG ) THEN
00923          IF( IPACK.EQ.1 ) THEN
00924 *
00925 *           'U' -- Upper triangular, not packed
00926 *
00927             DO 300 J = 1, M
00928                DO 290 I = J + 1, M
00929                   A( I, J ) = ZERO
00930   290          CONTINUE
00931   300       CONTINUE
00932 *
00933          ELSE IF( IPACK.EQ.2 ) THEN
00934 *
00935 *           'L' -- Lower triangular, not packed
00936 *
00937             DO 320 J = 2, M
00938                DO 310 I = 1, J - 1
00939                   A( I, J ) = ZERO
00940   310          CONTINUE
00941   320       CONTINUE
00942 *
00943          ELSE IF( IPACK.EQ.3 ) THEN
00944 *
00945 *           'C' -- Upper triangle packed Columnwise.
00946 *
00947             ICOL = 1
00948             IROW = 0
00949             DO 340 J = 1, M
00950                DO 330 I = 1, J
00951                   IROW = IROW + 1
00952                   IF( IROW.GT.LDA ) THEN
00953                      IROW = 1
00954                      ICOL = ICOL + 1
00955                   END IF
00956                   A( IROW, ICOL ) = A( I, J )
00957   330          CONTINUE
00958   340       CONTINUE
00959 *
00960          ELSE IF( IPACK.EQ.4 ) THEN
00961 *
00962 *           'R' -- Lower triangle packed Columnwise.
00963 *
00964             ICOL = 1
00965             IROW = 0
00966             DO 360 J = 1, M
00967                DO 350 I = J, M
00968                   IROW = IROW + 1
00969                   IF( IROW.GT.LDA ) THEN
00970                      IROW = 1
00971                      ICOL = ICOL + 1
00972                   END IF
00973                   A( IROW, ICOL ) = A( I, J )
00974   350          CONTINUE
00975   360       CONTINUE
00976 *
00977          ELSE IF( IPACK.GE.5 ) THEN
00978 *
00979 *           'B' -- The lower triangle is packed as a band matrix.
00980 *           'Q' -- The upper triangle is packed as a band matrix.
00981 *           'Z' -- The whole matrix is packed as a band matrix.
00982 *
00983             IF( IPACK.EQ.5 )
00984      $         UUB = 0
00985             IF( IPACK.EQ.6 )
00986      $         LLB = 0
00987 *
00988             DO 380 J = 1, UUB
00989                DO 370 I = MIN( J+LLB, M ), 1, -1
00990                   A( I-J+UUB+1, J ) = A( I, J )
00991   370          CONTINUE
00992   380       CONTINUE
00993 *
00994             DO 400 J = UUB + 2, N
00995                DO 390 I = J - UUB, MIN( J+LLB, M )
00996                   A( I-J+UUB+1, J ) = A( I, J )
00997   390          CONTINUE
00998   400       CONTINUE
00999          END IF
01000 *
01001 *        If packed, zero out extraneous elements.
01002 *
01003 *        Symmetric/Triangular Packed --
01004 *        zero out everything after A(IROW,ICOL)
01005 *
01006          IF( IPACK.EQ.3 .OR. IPACK.EQ.4 ) THEN
01007             DO 420 JC = ICOL, M
01008                DO 410 JR = IROW + 1, LDA
01009                   A( JR, JC ) = ZERO
01010   410          CONTINUE
01011                IROW = 0
01012   420       CONTINUE
01013 *
01014          ELSE IF( IPACK.GE.5 ) THEN
01015 *
01016 *           Packed Band --
01017 *              1st row is now in A( UUB+2-j, j), zero above it
01018 *              m-th row is now in A( M+UUB-j,j), zero below it
01019 *              last non-zero diagonal is now in A( UUB+LLB+1,j ),
01020 *                 zero below it, too.
01021 *
01022             IR1 = UUB + LLB + 2
01023             IR2 = UUB + M + 2
01024             DO 450 JC = 1, N
01025                DO 430 JR = 1, UUB + 1 - JC
01026                   A( JR, JC ) = ZERO
01027   430          CONTINUE
01028                DO 440 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
01029                   A( JR, JC ) = ZERO
01030   440          CONTINUE
01031   450       CONTINUE
01032          END IF
01033       END IF
01034 *
01035       RETURN
01036 *
01037 *     End of DLATMS
01038 *
01039       END