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