LAPACK: TESTING/MATGEN/clatms.f Source File

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