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