LAPACK: TESTING/MATGEN/clatmt.f Source File

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