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