LAPACK: TESTING/MATGEN/dlatmt.f Source File

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