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