|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE DLATMS( 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 DOUBLE PRECISION COND, DMAX 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER ISEED( 4 ) 00015 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLATMS generates random matrices with specified singular values 00022 * (or symmetric/hermitian with specified eigenvalues) 00023 * for testing LAPACK programs. 00024 * 00025 * DLATMS 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 DLATMS 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION 00147 * On entry, this is used as described under MODE above. 00148 * If used, it must be >= 1. Not modified. 00149 * 00150 * DMAX - DOUBLE PRECISION 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 DLATMS differ only in the PACK parameter, 00200 * they will generate mathematically equivalent matrices. 00201 * Not modified. 00202 * 00203 * A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 DLATM1 00247 * 2 => Cannot scale to DMAX (max. sing. value is 0) 00248 * 3 => Error return from DLAGGE or SLAGSY 00249 * 00250 * ===================================================================== 00251 * 00252 * .. Parameters .. 00253 DOUBLE PRECISION ZERO 00254 PARAMETER ( ZERO = 0.0D0 ) 00255 DOUBLE PRECISION ONE 00256 PARAMETER ( ONE = 1.0D0 ) 00257 DOUBLE PRECISION TWOPI 00258 PARAMETER ( TWOPI = 6.2831853071795864769252867663D+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 DOUBLE PRECISION ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP 00268 * .. 00269 * .. External Functions .. 00270 LOGICAL LSAME 00271 DOUBLE PRECISION DLARND 00272 EXTERNAL LSAME, DLARND 00273 * .. 00274 * .. External Subroutines .. 00275 EXTERNAL DCOPY, DLAGGE, DLAGSY, DLAROT, DLARTG, DLASET, 00276 $ DLATM1, DSCAL, XERBLA 00277 * .. 00278 * .. Intrinsic Functions .. 00279 INTRINSIC ABS, COS, DBLE, MAX, MIN, MOD, 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( DBLE( LLB+UUB ).LT.0.3D0*DBLE( 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( 'DLATMS', -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 DLATM1( 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 DSCAL( 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 DLASET( '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 DCOPY( 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 DCOPY( 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*DLARND( 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 DLAROT( .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 DLARTG( 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 DLAROT( .FALSE., ILTEMP, .TRUE., IL, C, -S, 00558 $ A( IROW-ISKEW*IC+IOFFST, IC ), 00559 $ ILDA, TEMP, EXTRA ) 00560 IF( ILTEMP ) THEN 00561 CALL DLARTG( 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 DLAROT( .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*DLARND( 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 DLAROT( .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 DLARTG( 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 DLAROT( .TRUE., ILTEMP, .TRUE., IL, C, -S, 00609 $ A( IR-ISKEW*ICOL+IOFFST, ICOL ), 00610 $ ILDA, TEMP, EXTRA ) 00611 IF( ILTEMP ) THEN 00612 CALL DLARTG( 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 DLAROT( .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*DLARND( 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 DLAROT( .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 DLARTG( 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 DLAROT( .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 DLARTG( 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 DLAROT( .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*DLARND( 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 DLAROT( .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 DLARTG( 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 DLAROT( .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 DLARTG( 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 DLAROT( .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 DCOPY( 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*DLARND( 1, ISEED ) 00765 C = COS( ANGLE ) 00766 S = SIN( ANGLE ) 00767 CALL DLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S, 00768 $ A( IROW-ISKEW*JC+IOFFG, JC ), ILDA, 00769 $ EXTRA, TEMP ) 00770 CALL DLAROT( .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 DLARTG( 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 DLAROT( .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 DLAROT( .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 DCOPY( 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*DLARND( 1, ISEED ) 00838 C = COS( ANGLE ) 00839 S = -SIN( ANGLE ) 00840 CALL DLAROT( .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 DLAROT( .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 DLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ), 00853 $ EXTRA, C, S, DUMMY ) 00854 TEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH ) 00855 CALL DLAROT( .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 DLAROT( .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 DLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK, 00908 $ IINFO ) 00909 ELSE 00910 * 00911 * Symmetric -- A = U D U' 00912 * 00913 CALL DLAGSY( 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 DLATMS 01040 * 01041 END
