Home
Browse Files
   1       SUBROUTINE DLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
   2      $                   KL, KU, PACK, A, LDA, WORK, INFO )
   3 *
   4 *  -- LAPACK test routine (version 3.1) --
   5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   6 *     June 2010
   7 *
   8 *     .. Scalar Arguments ..
   9       CHARACTER          DIST, PACK, SYM
  10       INTEGER            INFO, KL, KU, LDA, M, MODE, N
  11       DOUBLE PRECISION   COND, DMAX
  12 *     ..
  13 *     .. Array Arguments ..
  14       INTEGER            ISEED( 4 )
  15       DOUBLE PRECISION   A( LDA, * ), D( * ), WORK( * )
  16 *     ..
  17 *
  18 *  Purpose
  19 *  =======
  20 *
  21 *     DLATMS generates random matrices with specified singular values
  22 *     (or symmetric/hermitian with specified eigenvalues)
  23 *     for testing LAPACK programs.
  24 *
  25 *     DLATMS operates by applying the following sequence of
  26 *     operations:
  27 *
  28 *       Set the diagonal to D, where D may be input or
  29 *          computed according to MODE, COND, DMAX, and SYM
  30 *          as described below.
  31 *
  32 *       Generate a matrix with the appropriate band structure, by one
  33 *          of two methods:
  34 *
  35 *       Method A:
  36 *           Generate a dense M x N matrix by multiplying D on the left
  37 *               and the right by random unitary matrices, then:
  38 *
  39 *           Reduce the bandwidth according to KL and KU, using
  40 *           Householder transformations.
  41 *
  42 *       Method B:
  43 *           Convert the bandwidth-0 (i.e., diagonal) matrix to a
  44 *               bandwidth-1 matrix using Givens rotations, "chasing"
  45 *               out-of-band elements back, much as in QR; then
  46 *               convert the bandwidth-1 to a bandwidth-2 matrix, etc.
  47 *               Note that for reasonably small bandwidths (relative to
  48 *               M and N) this requires less storage, as a dense matrix
  49 *               is not generated.  Also, for symmetric matrices, only
  50 *               one triangle is generated.
  51 *
  52 *       Method A is chosen if the bandwidth is a large fraction of the
  53 *           order of the matrix, and LDA is at least M (so a dense
  54 *           matrix can be stored.)  Method B is chosen if the bandwidth
  55 *           is small (< 1/2 N for symmetric, < .3 N+M for
  56 *           non-symmetric), or LDA is less than M and not less than the
  57 *           bandwidth.
  58 *
  59 *       Pack the matrix if desired. Options specified by PACK are:
  60 *          no packing
  61 *          zero out upper half (if symmetric)
  62 *          zero out lower half (if symmetric)
  63 *          store the upper half columnwise (if symmetric or upper
  64 *                triangular)
  65 *          store the lower half columnwise (if symmetric or lower
  66 *                triangular)
  67 *          store the lower triangle in banded format (if symmetric
  68 *                or lower triangular)
  69 *          store the upper triangle in banded format (if symmetric
  70 *                or upper triangular)
  71 *          store the entire matrix in banded format
  72 *       If Method B is chosen, and band format is specified, then the
  73 *          matrix will be generated in the band format, so no repacking
  74 *          will be necessary.
  75 *
  76 *  Arguments
  77 *  =========
  78 *
  79 *  M        (input) INTEGER
  80 *           The number of rows of A. Not modified.
  81 *
  82 *  N        (input) INTEGER
  83 *           The number of columns of A. Not modified.
  84 *
  85 *  DIST     (input) CHARACTER*1
  86 *           On entry, DIST specifies the type of distribution to be used
  87 *           to generate the random eigen-/singular values.
  88 *           'U' => UNIFORM( 0, 1 )  ( 'U' for uniform )
  89 *           'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
  90 *           'N' => NORMAL( 0, 1 )   ( 'N' for normal )
  91 *           Not modified.
  92 *
  93 *  ISEED    (input/output) INTEGER array, dimension ( 4 )
  94 *           On entry ISEED specifies the seed of the random number
  95 *           generator. They should lie between 0 and 4095 inclusive,
  96 *           and ISEED(4) should be odd. The random number generator
  97 *           uses a linear congruential sequence limited to small
  98 *           integers, and so should produce machine independent
  99 *           random numbers. The values of ISEED are changed on
 100 *           exit, and can be used in the next call to DLATMS
 101 *           to continue the same random number sequence.
 102 *           Changed on exit.
 103 *
 104 *  SYM      (input) CHARACTER*1
 105 *           If SYM='S' or 'H', the generated matrix is symmetric, with
 106 *             eigenvalues specified by D, COND, MODE, and DMAX; they
 107 *             may be positive, negative, or zero.
 108 *           If SYM='P', the generated matrix is symmetric, with
 109 *             eigenvalues (= singular values) specified by D, COND,
 110 *             MODE, and DMAX; they will not be negative.
 111 *           If SYM='N', the generated matrix is nonsymmetric, with
 112 *             singular values specified by D, COND, MODE, and DMAX;
 113 *             they will not be negative.
 114 *           Not modified.
 115 *
 116 *  D        (input/output) DOUBLE PRECISION array, dimension ( MIN( M , N ) )
 117 *           This array is used to specify the singular values or
 118 *           eigenvalues of A (see SYM, above.)  If MODE=0, then D is
 119 *           assumed to contain the singular/eigenvalues, otherwise
 120 *           they will be computed according to MODE, COND, and DMAX,
 121 *           and placed in D.
 122 *           Modified if MODE is nonzero.
 123 *
 124 *  MODE     (input) INTEGER
 125 *           On entry this describes how the singular/eigenvalues are to
 126 *           be specified:
 127 *           MODE = 0 means use D as input
 128 *           MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
 129 *           MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
 130 *           MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
 131 *           MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
 132 *           MODE = 5 sets D to random numbers in the range
 133 *                    ( 1/COND , 1 ) such that their logarithms
 134 *                    are uniformly distributed.
 135 *           MODE = 6 set D to random numbers from same distribution
 136 *                    as the rest of the matrix.
 137 *           MODE < 0 has the same meaning as ABS(MODE), except that
 138 *              the order of the elements of D is reversed.
 139 *           Thus if MODE is positive, D has entries ranging from
 140 *              1 to 1/COND, if negative, from 1/COND to 1,
 141 *           If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
 142 *              the elements of D will also be multiplied by a random
 143 *              sign (i.e., +1 or -1.)
 144 *           Not modified.
 145 *
 146 *  COND     (input) DOUBLE PRECISION
 147 *           On entry, this is used as described under MODE above.
 148 *           If used, it must be >= 1. Not modified.
 149 *
 150 *  DMAX     (input) DOUBLE PRECISION
 151 *           If MODE is neither -6, 0 nor 6, the contents of D, as
 152 *           computed according to MODE and COND, will be scaled by
 153 *           DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
 154 *           singular value (which is to say the norm) will be abs(DMAX).
 155 *           Note that DMAX need not be positive: if DMAX is negative
 156 *           (or zero), D will be scaled by a negative number (or zero).
 157 *           Not modified.
 158 *
 159 *  KL       (input) INTEGER
 160 *           This specifies the lower bandwidth of the  matrix. For
 161 *           example, KL=0 implies upper triangular, KL=1 implies upper
 162 *           Hessenberg, and KL being at least M-1 means that the matrix
 163 *           has full lower bandwidth.  KL must equal KU if the matrix
 164 *           is symmetric.
 165 *           Not modified.
 166 *
 167 *  KU       (input) INTEGER
 168 *           This specifies the upper bandwidth of the  matrix. For
 169 *           example, KU=0 implies lower triangular, KU=1 implies lower
 170 *           Hessenberg, and KU being at least N-1 means that the matrix
 171 *           has full upper bandwidth.  KL must equal KU if the matrix
 172 *           is symmetric.
 173 *           Not modified.
 174 *
 175 *  PACK     (input) CHARACTER*1
 176 *           This specifies packing of matrix as follows:
 177 *           'N' => no packing
 178 *           'U' => zero out all subdiagonal entries (if symmetric)
 179 *           'L' => zero out all superdiagonal entries (if symmetric)
 180 *           'C' => store the upper triangle columnwise
 181 *                  (only if the matrix is symmetric or upper triangular)
 182 *           'R' => store the lower triangle columnwise
 183 *                  (only if the matrix is symmetric or lower triangular)
 184 *           'B' => store the lower triangle in band storage scheme
 185 *                  (only if matrix symmetric or lower triangular)
 186 *           'Q' => store the upper triangle in band storage scheme
 187 *                  (only if matrix symmetric or upper triangular)
 188 *           'Z' => store the entire matrix in band storage scheme
 189 *                      (pivoting can be provided for by using this
 190 *                      option to store A in the trailing rows of
 191 *                      the allocated storage)
 192 *
 193 *           Using these options, the various LAPACK packed and banded
 194 *           storage schemes can be obtained:
 195 *           GB               - use 'Z'
 196 *           PB, SB or TB     - use 'B' or 'Q'
 197 *           PP, SP or TP     - use 'C' or 'R'
 198 *
 199 *           If two calls to DLATMS differ only in the PACK parameter,
 200 *           they will generate mathematically equivalent matrices.
 201 *           Not modified.
 202 *
 203 *  A        (input/output) DOUBLE PRECISION array, dimension ( LDA, N )
 204 *           On exit A is the desired test matrix.  A is first generated
 205 *           in full (unpacked) form, and then packed, if so specified
 206 *           by PACK.  Thus, the first M elements of the first N
 207 *           columns will always be modified.  If PACK specifies a
 208 *           packed or banded storage scheme, all LDA elements of the
 209 *           first N columns will be modified; the elements of the
 210 *           array which do not correspond to elements of the generated
 211 *           matrix are set to zero.
 212 *           Modified.
 213 *
 214 *  LDA      (input) INTEGER
 215 *           LDA specifies the first dimension of A as declared in the
 216 *           calling program.  If PACK='N', 'U', 'L', 'C', or 'R', then
 217 *           LDA must be at least M.  If PACK='B' or 'Q', then LDA must
 218 *           be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
 219 *           If PACK='Z', LDA must be large enough to hold the packed
 220 *           array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
 221 *           Not modified.
 222 *
 223 *  WORK     (workspace) DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) )
 224 *           Workspace.
 225 *           Modified.
 226 *
 227 *  INFO     (output) INTEGER
 228 *           Error code.  On exit, INFO will be set to one of the
 229 *           following values:
 230 *             0 => normal return
 231 *            -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
 232 *            -2 => N negative
 233 *            -3 => DIST illegal string
 234 *            -5 => SYM illegal string
 235 *            -7 => MODE not in range -6 to 6
 236 *            -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
 237 *           -10 => KL negative
 238 *           -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
 239 *           -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
 240 *                  or PACK='C' or 'Q' and SYM='N' and KL is not zero;
 241 *                  or PACK='R' or 'B' and SYM='N' and KU is not zero;
 242 *                  or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
 243 *                  N.
 244 *           -14 => LDA is less than M, or PACK='Z' and LDA is less than
 245 *                  MIN(KU,N-1) + MIN(KL,M-1) + 1.
 246 *            1  => Error return from DLATM1
 247 *            2  => Cannot scale to DMAX (max. sing. value is 0)
 248 *            3  => Error return from DLAGGE or SLAGSY
 249 *
 250 *  =====================================================================
 251 *
 252 *     .. Parameters ..
 253       DOUBLE PRECISION   ZERO
 254       PARAMETER          ( ZERO = 0.0D0 )
 255       DOUBLE PRECISION   ONE
 256       PARAMETER          ( ONE = 1.0D0 )
 257       DOUBLE PRECISION   TWOPI
 258       PARAMETER          ( TWOPI = 6.2831853071795864769252867663D+0 )
 259 *     ..
 260 *     .. Local Scalars ..
 261       LOGICAL            GIVENS, ILEXTR, ILTEMP, TOPDWN
 262       INTEGER            I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
 263      $                   IOFFG, IOFFST, IPACK, IPACKG, IR, IR1, IR2,
 264      $                   IROW, IRSIGN, ISKEW, ISYM, ISYMPK, J, JC, JCH,
 265      $                   JKL, JKU, JR, K, LLB, MINLDA, MNMIN, MR, NC,
 266      $                   UUB
 267       DOUBLE PRECISION   ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
 268 *     ..
 269 *     .. External Functions ..
 270       LOGICAL            LSAME
 271       DOUBLE PRECISION   DLARND
 272       EXTERNAL           LSAME, DLARND
 273 *     ..
 274 *     .. External Subroutines ..
 275       EXTERNAL           DCOPY, DLAGGE, DLAGSY, DLAROT, DLARTG, DLASET,
 276      $                   DLATM1, DSCAL, XERBLA
 277 *     ..
 278 *     .. Intrinsic Functions ..
 279       INTRINSIC          ABS, COS, DBLE, MAX, MIN, MOD, SIN
 280 *     ..
 281 *     .. Executable Statements ..
 282 *
 283 *     1)      Decode and Test the input parameters.
 284 *             Initialize flags & seed.
 285 *
 286       INFO = 0
 287 *
 288 *     Quick return if possible
 289 *
 290       IF( M.EQ.0 .OR. N.EQ.0 )
 291      $   RETURN
 292 *
 293 *     Decode DIST
 294 *
 295       IF( LSAME( DIST, 'U' ) ) THEN
 296          IDIST = 1
 297       ELSE IF( LSAME( DIST, 'S' ) ) THEN
 298          IDIST = 2
 299       ELSE IF( LSAME( DIST, 'N' ) ) THEN
 300          IDIST = 3
 301       ELSE
 302          IDIST = -1
 303       END IF
 304 *
 305 *     Decode SYM
 306 *
 307       IF( LSAME( SYM, 'N' ) ) THEN
 308          ISYM = 1
 309          IRSIGN = 0
 310       ELSE IF( LSAME( SYM, 'P' ) ) THEN
 311          ISYM = 2
 312          IRSIGN = 0
 313       ELSE IF( LSAME( SYM, 'S' ) ) THEN
 314          ISYM = 2
 315          IRSIGN = 1
 316       ELSE IF( LSAME( SYM, 'H' ) ) THEN
 317          ISYM = 2
 318          IRSIGN = 1
 319       ELSE
 320          ISYM = -1
 321       END IF
 322 *
 323 *     Decode PACK
 324 *
 325       ISYMPK = 0
 326       IF( LSAME( PACK, 'N' ) ) THEN
 327          IPACK = 0
 328       ELSE IF( LSAME( PACK, 'U' ) ) THEN
 329          IPACK = 1
 330          ISYMPK = 1
 331       ELSE IF( LSAME( PACK, 'L' ) ) THEN
 332          IPACK = 2
 333          ISYMPK = 1
 334       ELSE IF( LSAME( PACK, 'C' ) ) THEN
 335          IPACK = 3
 336          ISYMPK = 2
 337       ELSE IF( LSAME( PACK, 'R' ) ) THEN
 338          IPACK = 4
 339          ISYMPK = 3
 340       ELSE IF( LSAME( PACK, 'B' ) ) THEN
 341          IPACK = 5
 342          ISYMPK = 3
 343       ELSE IF( LSAME( PACK, 'Q' ) ) THEN
 344          IPACK = 6
 345          ISYMPK = 2
 346       ELSE IF( LSAME( PACK, 'Z' ) ) THEN
 347          IPACK = 7
 348       ELSE
 349          IPACK = -1
 350       END IF
 351 *
 352 *     Set certain internal parameters
 353 *
 354       MNMIN = MIN( M, N )
 355       LLB = MIN( KL, M-1 )
 356       UUB = MIN( KU, N-1 )
 357       MR = MIN( M, N+LLB )
 358       NC = MIN( N, M+UUB )
 359 *
 360       IF( IPACK.EQ.5 .OR. IPACK.EQ.6 ) THEN
 361          MINLDA = UUB + 1
 362       ELSE IF( IPACK.EQ.7 ) THEN
 363          MINLDA = LLB + UUB + 1
 364       ELSE
 365          MINLDA = M
 366       END IF
 367 *
 368 *     Use Givens rotation method if bandwidth small enough,
 369 *     or if LDA is too small to store the matrix unpacked.
 370 *
 371       GIVENS = .FALSE.
 372       IF( ISYM.EQ.1 ) THEN
 373          IF( DBLE( LLB+UUB ).LT.0.3D0*DBLE( MAX( 1, MR+NC ) ) )
 374      $      GIVENS = .TRUE.
 375       ELSE
 376          IF( 2*LLB.LT.M )
 377      $      GIVENS = .TRUE.
 378       END IF
 379       IF( LDA.LT.M .AND. LDA.GE.MINLDA )
 380      $   GIVENS = .TRUE.
 381 *
 382 *     Set INFO if an error
 383 *
 384       IF( M.LT.0 ) THEN
 385          INFO = -1
 386       ELSE IF( M.NE.N .AND. ISYM.NE.1 ) THEN
 387          INFO = -1
 388       ELSE IF( N.LT.0 ) THEN
 389          INFO = -2
 390       ELSE IF( IDIST.EQ.-1 ) THEN
 391          INFO = -3
 392       ELSE IF( ISYM.EQ.-1 ) THEN
 393          INFO = -5
 394       ELSE IF( ABS( MODE ).GT.6 ) THEN
 395          INFO = -7
 396       ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
 397      $          THEN
 398          INFO = -8
 399       ELSE IF( KL.LT.0 ) THEN
 400          INFO = -10
 401       ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN
 402          INFO = -11
 403       ELSE IF( IPACK.EQ.-1 .OR. ( ISYMPK.EQ.1 .AND. ISYM.EQ.1 ) .OR.
 404      $         ( ISYMPK.EQ.2 .AND. ISYM.EQ.1 .AND. KL.GT.0 ) .OR.
 405      $         ( ISYMPK.EQ.3 .AND. ISYM.EQ.1 .AND. KU.GT.0 ) .OR.
 406      $         ( ISYMPK.NE.0 .AND. M.NE.N ) ) THEN
 407          INFO = -12
 408       ELSE IF( LDA.LT.MAX( 1, MINLDA ) ) THEN
 409          INFO = -14
 410       END IF
 411 *
 412       IF( INFO.NE.0 ) THEN
 413          CALL XERBLA( 'DLATMS', -INFO )
 414          RETURN
 415       END IF
 416 *
 417 *     Initialize random number generator
 418 *
 419       DO 10 I = 1, 4
 420          ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
 421    10 CONTINUE
 422 *
 423       IF( MOD( ISEED( 4 ), 2 ).NE.1 )
 424      $   ISEED( 4 ) = ISEED( 4 ) + 1
 425 *
 426 *     2)      Set up D  if indicated.
 427 *
 428 *             Compute D according to COND and MODE
 429 *
 430       CALL DLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, IINFO )
 431       IF( IINFO.NE.0 ) THEN
 432          INFO = 1
 433          RETURN
 434       END IF
 435 *
 436 *     Choose Top-Down if D is (apparently) increasing,
 437 *     Bottom-Up if D is (apparently) decreasing.
 438 *
 439       IF( ABS( D( 1 ) ).LE.ABS( D( MNMIN ) ) ) THEN
 440          TOPDWN = .TRUE.
 441       ELSE
 442          TOPDWN = .FALSE.
 443       END IF
 444 *
 445       IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
 446 *
 447 *        Scale by DMAX
 448 *
 449          TEMP = ABS( D( 1 ) )
 450          DO 20 I = 2, MNMIN
 451             TEMP = MAX( TEMP, ABS( D( I ) ) )
 452    20    CONTINUE
 453 *
 454          IF( TEMP.GT.ZERO ) THEN
 455             ALPHA = DMAX / TEMP
 456          ELSE
 457             INFO = 2
 458             RETURN
 459          END IF
 460 *
 461          CALL DSCAL( MNMIN, ALPHA, D, 1 )
 462 *
 463       END IF
 464 *
 465 *     3)      Generate Banded Matrix using Givens rotations.
 466 *             Also the special case of UUB=LLB=0
 467 *
 468 *               Compute Addressing constants to cover all
 469 *               storage formats.  Whether GE, SY, GB, or SB,
 470 *               upper or lower triangle or both,
 471 *               the (i,j)-th element is in
 472 *               A( i - ISKEW*j + IOFFST, j )
 473 *
 474       IF( IPACK.GT.4 ) THEN
 475          ILDA = LDA - 1
 476          ISKEW = 1
 477          IF( IPACK.GT.5 ) THEN
 478             IOFFST = UUB + 1
 479          ELSE
 480             IOFFST = 1
 481          END IF
 482       ELSE
 483          ILDA = LDA
 484          ISKEW = 0
 485          IOFFST = 0
 486       END IF
 487 *
 488 *     IPACKG is the format that the matrix is generated in. If this is
 489 *     different from IPACK, then the matrix must be repacked at the
 490 *     end.  It also signals how to compute the norm, for scaling.
 491 *
 492       IPACKG = 0
 493       CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
 494 *
 495 *     Diagonal Matrix -- We are done, unless it
 496 *     is to be stored SP/PP/TP (PACK='R' or 'C')
 497 *
 498       IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN
 499          CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
 500          IF( IPACK.LE.2 .OR. IPACK.GE.5 )
 501      $      IPACKG = IPACK
 502 *
 503       ELSE IF( GIVENS ) THEN
 504 *
 505 *        Check whether to use Givens rotations,
 506 *        Householder transformations, or nothing.
 507 *
 508          IF( ISYM.EQ.1 ) THEN
 509 *
 510 *           Non-symmetric -- A = U D V
 511 *
 512             IF( IPACK.GT.4 ) THEN
 513                IPACKG = IPACK
 514             ELSE
 515                IPACKG = 0
 516             END IF
 517 *
 518             CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFST, 1 ), ILDA+1 )
 519 *
 520             IF( TOPDWN ) THEN
 521                JKL = 0
 522                DO 50 JKU = 1, UUB
 523 *
 524 *                 Transform from bandwidth JKL, JKU-1 to JKL, JKU
 525 *
 526 *                 Last row actually rotated is M
 527 *                 Last column actually rotated is MIN( M+JKU, N )
 528 *
 529                   DO 40 JR = 1, MIN( M+JKU, N ) + JKL - 1
 530                      EXTRA = ZERO
 531                      ANGLE = TWOPI*DLARND( 1, ISEED )
 532                      C = COS( ANGLE )
 533                      S = SIN( ANGLE )
 534                      ICOL = MAX( 1, JR-JKL )
 535                      IF( JR.LT.M ) THEN
 536                         IL = MIN( N, JR+JKU ) + 1 - ICOL
 537                         CALL DLAROT( .TRUE., JR.GT.JKL, .FALSE., IL, C,
 538      $                               S, A( JR-ISKEW*ICOL+IOFFST, ICOL ),
 539      $                               ILDA, EXTRA, DUMMY )
 540                      END IF
 541 *
 542 *                    Chase "EXTRA" back up
 543 *
 544                      IR = JR
 545                      IC = ICOL
 546                      DO 30 JCH = JR - JKL, 1, -JKL - JKU
 547                         IF( IR.LT.M ) THEN
 548                            CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
 549      $                                  IC+1 ), EXTRA, C, S, DUMMY )
 550                         END IF
 551                         IROW = MAX( 1, JCH-JKU )
 552                         IL = IR + 2 - IROW
 553                         TEMP = ZERO
 554                         ILTEMP = JCH.GT.JKU
 555                         CALL DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S,
 556      $                               A( IROW-ISKEW*IC+IOFFST, IC ),
 557      $                               ILDA, TEMP, EXTRA )
 558                         IF( ILTEMP ) THEN
 559                            CALL DLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST,
 560      $                                  IC+1 ), TEMP, C, S, DUMMY )
 561                            ICOL = MAX( 1, JCH-JKU-JKL )
 562                            IL = IC + 2 - ICOL
 563                            EXTRA = ZERO
 564                            CALL DLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE.,
 565      $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
 566      $                                  IOFFST, ICOL ), ILDA, EXTRA,
 567      $                                  TEMP )
 568                            IC = ICOL
 569                            IR = IROW
 570                         END IF
 571    30                CONTINUE
 572    40             CONTINUE
 573    50          CONTINUE
 574 *
 575                JKU = UUB
 576                DO 80 JKL = 1, LLB
 577 *
 578 *                 Transform from bandwidth JKL-1, JKU to JKL, JKU
 579 *
 580                   DO 70 JC = 1, MIN( N+JKL, M ) + JKU - 1
 581                      EXTRA = ZERO
 582                      ANGLE = TWOPI*DLARND( 1, ISEED )
 583                      C = COS( ANGLE )
 584                      S = SIN( ANGLE )
 585                      IROW = MAX( 1, JC-JKU )
 586                      IF( JC.LT.N ) THEN
 587                         IL = MIN( M, JC+JKL ) + 1 - IROW
 588                         CALL DLAROT( .FALSE., JC.GT.JKU, .FALSE., IL, C,
 589      $                               S, A( IROW-ISKEW*JC+IOFFST, JC ),
 590      $                               ILDA, EXTRA, DUMMY )
 591                      END IF
 592 *
 593 *                    Chase "EXTRA" back up
 594 *
 595                      IC = JC
 596                      IR = IROW
 597                      DO 60 JCH = JC - JKU, 1, -JKL - JKU
 598                         IF( IC.LT.N ) THEN
 599                            CALL DLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST,
 600      $                                  IC+1 ), EXTRA, C, S, DUMMY )
 601                         END IF
 602                         ICOL = MAX( 1, JCH-JKL )
 603                         IL = IC + 2 - ICOL
 604                         TEMP = ZERO
 605                         ILTEMP = JCH.GT.JKL
 606                         CALL DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S,
 607      $                               A( IR-ISKEW*ICOL+IOFFST, ICOL ),
 608      $                               ILDA, TEMP, EXTRA )
 609                         IF( ILTEMP ) THEN
 610                            CALL DLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST,
 611      $                                  ICOL+1 ), TEMP, C, S, DUMMY )
 612                            IROW = MAX( 1, JCH-JKL-JKU )
 613                            IL = IR + 2 - IROW
 614                            EXTRA = ZERO
 615                            CALL DLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE.,
 616      $                                  IL, C, -S, A( IROW-ISKEW*ICOL+
 617      $                                  IOFFST, ICOL ), ILDA, EXTRA,
 618      $                                  TEMP )
 619                            IC = ICOL
 620                            IR = IROW
 621                         END IF
 622    60                CONTINUE
 623    70             CONTINUE
 624    80          CONTINUE
 625 *
 626             ELSE
 627 *
 628 *              Bottom-Up -- Start at the bottom right.
 629 *
 630                JKL = 0
 631                DO 110 JKU = 1, UUB
 632 *
 633 *                 Transform from bandwidth JKL, JKU-1 to JKL, JKU
 634 *
 635 *                 First row actually rotated is M
 636 *                 First column actually rotated is MIN( M+JKU, N )
 637 *
 638                   IENDCH = MIN( M, N+JKL ) - 1
 639                   DO 100 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1
 640                      EXTRA = ZERO
 641                      ANGLE = TWOPI*DLARND( 1, ISEED )
 642                      C = COS( ANGLE )
 643                      S = SIN( ANGLE )
 644                      IROW = MAX( 1, JC-JKU+1 )
 645                      IF( JC.GT.0 ) THEN
 646                         IL = MIN( M, JC+JKL+1 ) + 1 - IROW
 647                         CALL DLAROT( .FALSE., .FALSE., JC+JKL.LT.M, IL,
 648      $                               C, S, A( IROW-ISKEW*JC+IOFFST,
 649      $                               JC ), ILDA, DUMMY, EXTRA )
 650                      END IF
 651 *
 652 *                    Chase "EXTRA" back down
 653 *
 654                      IC = JC
 655                      DO 90 JCH = JC + JKL, IENDCH, JKL + JKU
 656                         ILEXTR = IC.GT.0
 657                         IF( ILEXTR ) THEN
 658                            CALL DLARTG( A( JCH-ISKEW*IC+IOFFST, IC ),
 659      $                                  EXTRA, C, S, DUMMY )
 660                         END IF
 661                         IC = MAX( 1, IC )
 662                         ICOL = MIN( N-1, JCH+JKU )
 663                         ILTEMP = JCH + JKU.LT.N
 664                         TEMP = ZERO
 665                         CALL DLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC,
 666      $                               C, S, A( JCH-ISKEW*IC+IOFFST, IC ),
 667      $                               ILDA, EXTRA, TEMP )
 668                         IF( ILTEMP ) THEN
 669                            CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFST,
 670      $                                  ICOL ), TEMP, C, S, DUMMY )
 671                            IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
 672                            EXTRA = ZERO
 673                            CALL DLAROT( .FALSE., .TRUE.,
 674      $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
 675      $                                  A( JCH-ISKEW*ICOL+IOFFST,
 676      $                                  ICOL ), ILDA, TEMP, EXTRA )
 677                            IC = ICOL
 678                         END IF
 679    90                CONTINUE
 680   100             CONTINUE
 681   110          CONTINUE
 682 *
 683                JKU = UUB
 684                DO 140 JKL = 1, LLB
 685 *
 686 *                 Transform from bandwidth JKL-1, JKU to JKL, JKU
 687 *
 688 *                 First row actually rotated is MIN( N+JKL, M )
 689 *                 First column actually rotated is N
 690 *
 691                   IENDCH = MIN( N, M+JKU ) - 1
 692                   DO 130 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1
 693                      EXTRA = ZERO
 694                      ANGLE = TWOPI*DLARND( 1, ISEED )
 695                      C = COS( ANGLE )
 696                      S = SIN( ANGLE )
 697                      ICOL = MAX( 1, JR-JKL+1 )
 698                      IF( JR.GT.0 ) THEN
 699                         IL = MIN( N, JR+JKU+1 ) + 1 - ICOL
 700                         CALL DLAROT( .TRUE., .FALSE., JR+JKU.LT.N, IL,
 701      $                               C, S, A( JR-ISKEW*ICOL+IOFFST,
 702      $                               ICOL ), ILDA, DUMMY, EXTRA )
 703                      END IF
 704 *
 705 *                    Chase "EXTRA" back down
 706 *
 707                      IR = JR
 708                      DO 120 JCH = JR + JKU, IENDCH, JKL + JKU
 709                         ILEXTR = IR.GT.0
 710                         IF( ILEXTR ) THEN
 711                            CALL DLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ),
 712      $                                  EXTRA, C, S, DUMMY )
 713                         END IF
 714                         IR = MAX( 1, IR )
 715                         IROW = MIN( M-1, JCH+JKL )
 716                         ILTEMP = JCH + JKL.LT.M
 717                         TEMP = ZERO
 718                         CALL DLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR,
 719      $                               C, S, A( IR-ISKEW*JCH+IOFFST,
 720      $                               JCH ), ILDA, EXTRA, TEMP )
 721                         IF( ILTEMP ) THEN
 722                            CALL DLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ),
 723      $                                  TEMP, C, S, DUMMY )
 724                            IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH
 725                            EXTRA = ZERO
 726                            CALL DLAROT( .TRUE., .TRUE.,
 727      $                                  JCH+JKL+JKU.LE.IENDCH, IL, C, S,
 728      $                                  A( IROW-ISKEW*JCH+IOFFST, JCH ),
 729      $                                  ILDA, TEMP, EXTRA )
 730                            IR = IROW
 731                         END IF
 732   120                CONTINUE
 733   130             CONTINUE
 734   140          CONTINUE
 735             END IF
 736 *
 737          ELSE
 738 *
 739 *           Symmetric -- A = U D U'
 740 *
 741             IPACKG = IPACK
 742             IOFFG = IOFFST
 743 *
 744             IF( TOPDWN ) THEN
 745 *
 746 *              Top-Down -- Generate Upper triangle only
 747 *
 748                IF( IPACK.GE.5 ) THEN
 749                   IPACKG = 6
 750                   IOFFG = UUB + 1
 751                ELSE
 752                   IPACKG = 1
 753                END IF
 754                CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
 755 *
 756                DO 170 K = 1, UUB
 757                   DO 160 JC = 1, N - 1
 758                      IROW = MAX( 1, JC-K )
 759                      IL = MIN( JC+1, K+2 )
 760                      EXTRA = ZERO
 761                      TEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 )
 762                      ANGLE = TWOPI*DLARND( 1, ISEED )
 763                      C = COS( ANGLE )
 764                      S = SIN( ANGLE )
 765                      CALL DLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S,
 766      $                            A( IROW-ISKEW*JC+IOFFG, JC ), ILDA,
 767      $                            EXTRA, TEMP )
 768                      CALL DLAROT( .TRUE., .TRUE., .FALSE.,
 769      $                            MIN( K, N-JC )+1, C, S,
 770      $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
 771      $                            TEMP, DUMMY )
 772 *
 773 *                    Chase EXTRA back up the matrix
 774 *
 775                      ICOL = JC
 776                      DO 150 JCH = JC - K, 1, -K
 777                         CALL DLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG,
 778      $                               ICOL+1 ), EXTRA, C, S, DUMMY )
 779                         TEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 )
 780                         CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, -S,
 781      $                               A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
 782      $                               ILDA, TEMP, EXTRA )
 783                         IROW = MAX( 1, JCH-K )
 784                         IL = MIN( JCH+1, K+2 )
 785                         EXTRA = ZERO
 786                         CALL DLAROT( .FALSE., JCH.GT.K, .TRUE., IL, C,
 787      $                               -S, A( IROW-ISKEW*JCH+IOFFG, JCH ),
 788      $                               ILDA, EXTRA, TEMP )
 789                         ICOL = JCH
 790   150                CONTINUE
 791   160             CONTINUE
 792   170          CONTINUE
 793 *
 794 *              If we need lower triangle, copy from upper. Note that
 795 *              the order of copying is chosen to work for 'q' -> 'b'
 796 *
 797                IF( IPACK.NE.IPACKG .AND. IPACK.NE.3 ) THEN
 798                   DO 190 JC = 1, N
 799                      IROW = IOFFST - ISKEW*JC
 800                      DO 180 JR = JC, MIN( N, JC+UUB )
 801                         A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
 802   180                CONTINUE
 803   190             CONTINUE
 804                   IF( IPACK.EQ.5 ) THEN
 805                      DO 210 JC = N - UUB + 1, N
 806                         DO 200 JR = N + 2 - JC, UUB + 1
 807                            A( JR, JC ) = ZERO
 808   200                   CONTINUE
 809   210                CONTINUE
 810                   END IF
 811                   IF( IPACKG.EQ.6 ) THEN
 812                      IPACKG = IPACK
 813                   ELSE
 814                      IPACKG = 0
 815                   END IF
 816                END IF
 817             ELSE
 818 *
 819 *              Bottom-Up -- Generate Lower triangle only
 820 *
 821                IF( IPACK.GE.5 ) THEN
 822                   IPACKG = 5
 823                   IF( IPACK.EQ.6 )
 824      $               IOFFG = 1
 825                ELSE
 826                   IPACKG = 2
 827                END IF
 828                CALL DCOPY( MNMIN, D, 1, A( 1-ISKEW+IOFFG, 1 ), ILDA+1 )
 829 *
 830                DO 240 K = 1, UUB
 831                   DO 230 JC = N - 1, 1, -1
 832                      IL = MIN( N+1-JC, K+2 )
 833                      EXTRA = ZERO
 834                      TEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC )
 835                      ANGLE = TWOPI*DLARND( 1, ISEED )
 836                      C = COS( ANGLE )
 837                      S = -SIN( ANGLE )
 838                      CALL DLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S,
 839      $                            A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA,
 840      $                            TEMP, EXTRA )
 841                      ICOL = MAX( 1, JC-K+1 )
 842                      CALL DLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, C,
 843      $                            S, A( JC-ISKEW*ICOL+IOFFG, ICOL ),
 844      $                            ILDA, DUMMY, TEMP )
 845 *
 846 *                    Chase EXTRA back down the matrix
 847 *
 848                      ICOL = JC
 849                      DO 220 JCH = JC + K, N - 1, K
 850                         CALL DLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
 851      $                               EXTRA, C, S, DUMMY )
 852                         TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH )
 853                         CALL DLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S,
 854      $                               A( JCH-ISKEW*ICOL+IOFFG, ICOL ),
 855      $                               ILDA, EXTRA, TEMP )
 856                         IL = MIN( N+1-JCH, K+2 )
 857                         EXTRA = ZERO
 858                         CALL DLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, C,
 859      $                               S, A( ( 1-ISKEW )*JCH+IOFFG, JCH ),
 860      $                               ILDA, TEMP, EXTRA )
 861                         ICOL = JCH
 862   220                CONTINUE
 863   230             CONTINUE
 864   240          CONTINUE
 865 *
 866 *              If we need upper triangle, copy from lower. Note that
 867 *              the order of copying is chosen to work for 'b' -> 'q'
 868 *
 869                IF( IPACK.NE.IPACKG .AND. IPACK.NE.4 ) THEN
 870                   DO 260 JC = N, 1, -1
 871                      IROW = IOFFST - ISKEW*JC
 872                      DO 250 JR = JC, MAX( 1, JC-UUB ), -1
 873                         A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR )
 874   250                CONTINUE
 875   260             CONTINUE
 876                   IF( IPACK.EQ.6 ) THEN
 877                      DO 280 JC = 1, UUB
 878                         DO 270 JR = 1, UUB + 1 - JC
 879                            A( JR, JC ) = ZERO
 880   270                   CONTINUE
 881   280                CONTINUE
 882                   END IF
 883                   IF( IPACKG.EQ.5 ) THEN
 884                      IPACKG = IPACK
 885                   ELSE
 886                      IPACKG = 0
 887                   END IF
 888                END IF
 889             END IF
 890          END IF
 891 *
 892       ELSE
 893 *
 894 *        4)      Generate Banded Matrix by first
 895 *                Rotating by random Unitary matrices,
 896 *                then reducing the bandwidth using Householder
 897 *                transformations.
 898 *
 899 *                Note: we should get here only if LDA .ge. N
 900 *
 901          IF( ISYM.EQ.1 ) THEN
 902 *
 903 *           Non-symmetric -- A = U D V
 904 *
 905             CALL DLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK,
 906      $                   IINFO )
 907          ELSE
 908 *
 909 *           Symmetric -- A = U D U'
 910 *
 911             CALL DLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO )
 912 *
 913          END IF
 914          IF( IINFO.NE.0 ) THEN
 915             INFO = 3
 916             RETURN
 917          END IF
 918       END IF
 919 *
 920 *     5)      Pack the matrix
 921 *
 922       IF( IPACK.NE.IPACKG ) THEN
 923          IF( IPACK.EQ.1 ) THEN
 924 *
 925 *           'U' -- Upper triangular, not packed
 926 *
 927             DO 300 J = 1, M
 928                DO 290 I = J + 1, M
 929                   A( I, J ) = ZERO
 930   290          CONTINUE
 931   300       CONTINUE
 932 *
 933          ELSE IF( IPACK.EQ.2 ) THEN
 934 *
 935 *           'L' -- Lower triangular, not packed
 936 *
 937             DO 320 J = 2, M
 938                DO 310 I = 1, J - 1
 939                   A( I, J ) = ZERO
 940   310          CONTINUE
 941   320       CONTINUE
 942 *
 943          ELSE IF( IPACK.EQ.3 ) THEN
 944 *
 945 *           'C' -- Upper triangle packed Columnwise.
 946 *
 947             ICOL = 1
 948             IROW = 0
 949             DO 340 J = 1, M
 950                DO 330 I = 1, J
 951                   IROW = IROW + 1
 952                   IF( IROW.GT.LDA ) THEN
 953                      IROW = 1
 954                      ICOL = ICOL + 1
 955                   END IF
 956                   A( IROW, ICOL ) = A( I, J )
 957   330          CONTINUE
 958   340       CONTINUE
 959 *
 960          ELSE IF( IPACK.EQ.4 ) THEN
 961 *
 962 *           'R' -- Lower triangle packed Columnwise.
 963 *
 964             ICOL = 1
 965             IROW = 0
 966             DO 360 J = 1, M
 967                DO 350 I = J, M
 968                   IROW = IROW + 1
 969                   IF( IROW.GT.LDA ) THEN
 970                      IROW = 1
 971                      ICOL = ICOL + 1
 972                   END IF
 973                   A( IROW, ICOL ) = A( I, J )
 974   350          CONTINUE
 975   360       CONTINUE
 976 *
 977          ELSE IF( IPACK.GE.5 ) THEN
 978 *
 979 *           'B' -- The lower triangle is packed as a band matrix.
 980 *           'Q' -- The upper triangle is packed as a band matrix.
 981 *           'Z' -- The whole matrix is packed as a band matrix.
 982 *
 983             IF( IPACK.EQ.5 )
 984      $         UUB = 0
 985             IF( IPACK.EQ.6 )
 986      $         LLB = 0
 987 *
 988             DO 380 J = 1, UUB
 989                DO 370 I = MIN( J+LLB, M ), 1, -1
 990                   A( I-J+UUB+1, J ) = A( I, J )
 991   370          CONTINUE
 992   380       CONTINUE
 993 *
 994             DO 400 J = UUB + 2, N
 995                DO 390 I = J - UUB, MIN( J+LLB, M )
 996                   A( I-J+UUB+1, J ) = A( I, J )
 997   390          CONTINUE
 998   400       CONTINUE
 999          END IF
1000 *
1001 *        If packed, zero out extraneous elements.
1002 *
1003 *        Symmetric/Triangular Packed --
1004 *        zero out everything after A(IROW,ICOL)
1005 *
1006          IF( IPACK.EQ.3 .OR. IPACK.EQ.4 ) THEN
1007             DO 420 JC = ICOL, M
1008                DO 410 JR = IROW + 1, LDA
1009                   A( JR, JC ) = ZERO
1010   410          CONTINUE
1011                IROW = 0
1012   420       CONTINUE
1013 *
1014          ELSE IF( IPACK.GE.5 ) THEN
1015 *
1016 *           Packed Band --
1017 *              1st row is now in A( UUB+2-j, j), zero above it
1018 *              m-th row is now in A( M+UUB-j,j), zero below it
1019 *              last non-zero diagonal is now in A( UUB+LLB+1,j ),
1020 *                 zero below it, too.
1021 *
1022             IR1 = UUB + LLB + 2
1023             IR2 = UUB + M + 2
1024             DO 450 JC = 1, N
1025                DO 430 JR = 1, UUB + 1 - JC
1026                   A( JR, JC ) = ZERO
1027   430          CONTINUE
1028                DO 440 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA
1029                   A( JR, JC ) = ZERO
1030   440          CONTINUE
1031   450       CONTINUE
1032          END IF
1033       END IF
1034 *
1035       RETURN
1036 *
1037 *     End of DLATMS
1038 *
1039       END
Home