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