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