|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE CLATMS( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, 00002 $ KL, KU, PACK, A, LDA, WORK, INFO ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER DIST, PACK, SYM 00010 INTEGER INFO, KL, KU, LDA, M, MODE, N 00011 REAL COND, DMAX 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER ISEED( 4 ) 00015 REAL D( * ) 00016 COMPLEX A( LDA, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CLATMS generates random matrices with specified singular values 00023 * (or hermitian with specified eigenvalues) 00024 * for testing LAPACK programs. 00025 * 00026 * CLATMS 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 - INTEGER 00081 * The number of rows of A. Not modified. 00082 * 00083 * N - 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 - 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 - 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 CLATMS 00104 * to continue the same random number sequence. 00105 * Changed on exit. 00106 * 00107 * SYM - 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 - REAL 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 - 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 - REAL 00153 * On entry, this is used as described under MODE above. 00154 * If used, it must be >= 1. Not modified. 00155 * 00156 * DMAX - REAL 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 - 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 - 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 - 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 CLATMS differ only in the PACK parameter, 00210 * they will generate mathematically equivalent matrices. 00211 * Not modified. 00212 * 00213 * A - COMPLEX 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 - 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 - COMPLEX array, dimension ( 3*MAX( N, M ) ) 00234 * Workspace. 00235 * Modified. 00236 * 00237 * INFO - 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 SLATM1 00258 * 2 => Cannot scale to DMAX (max. sing. value is 0) 00259 * 3 => Error return from CLAGGE, CLAGHE or CLAGSY 00260 * 00261 * ===================================================================== 00262 * 00263 * .. Parameters .. 00264 REAL ZERO 00265 PARAMETER ( ZERO = 0.0E+0 ) 00266 REAL ONE 00267 PARAMETER ( ONE = 1.0E+0 ) 00268 COMPLEX CZERO 00269 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00270 REAL TWOPI 00271 PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 ) 00272 * .. 00273 * .. Local Scalars .. 00274 LOGICAL CSYM, GIVENS, ILEXTR, ILTEMP, TOPDWN 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 REAL ALPHA, ANGLE, REALC, TEMP 00281 COMPLEX C, CT, CTEMP, DUMMY, EXTRA, S, ST 00282 * .. 00283 * .. External Functions .. 00284 LOGICAL LSAME 00285 REAL SLARND 00286 COMPLEX CLARND 00287 EXTERNAL LSAME, SLARND, CLARND 00288 * .. 00289 * .. External Subroutines .. 00290 EXTERNAL CLAGGE, CLAGHE, CLAGSY, CLAROT, CLARTG, CLASET, 00291 $ SLATM1, SSCAL, XERBLA 00292 * .. 00293 * .. Intrinsic Functions .. 00294 INTRINSIC ABS, CMPLX, CONJG, COS, MAX, MIN, MOD, REAL, 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 CSYM = .FALSE. 00327 ELSE IF( LSAME( SYM, 'P' ) ) THEN 00328 ISYM = 2 00329 IRSIGN = 0 00330 CSYM = .FALSE. 00331 ELSE IF( LSAME( SYM, 'S' ) ) THEN 00332 ISYM = 2 00333 IRSIGN = 0 00334 CSYM = .TRUE. 00335 ELSE IF( LSAME( SYM, 'H' ) ) THEN 00336 ISYM = 2 00337 IRSIGN = 1 00338 CSYM = .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 IROW = 1 00380 ICOL = 1 00381 CSYM = .FALSE. 00382 * 00383 IF( IPACK.EQ.5 .OR. IPACK.EQ.6 ) THEN 00384 MINLDA = UUB + 1 00385 ELSE IF( IPACK.EQ.7 ) THEN 00386 MINLDA = LLB + UUB + 1 00387 ELSE 00388 MINLDA = M 00389 END IF 00390 * 00391 * Use Givens rotation method if bandwidth small enough, 00392 * or if LDA is too small to store the matrix unpacked. 00393 * 00394 GIVENS = .FALSE. 00395 IF( ISYM.EQ.1 ) THEN 00396 IF( REAL( LLB+UUB ).LT.0.3*REAL( MAX( 1, MR+NC ) ) ) 00397 $ GIVENS = .TRUE. 00398 ELSE 00399 IF( 2*LLB.LT.M ) 00400 $ GIVENS = .TRUE. 00401 END IF 00402 IF( LDA.LT.M .AND. LDA.GE.MINLDA ) 00403 $ GIVENS = .TRUE. 00404 * 00405 * Set INFO if an error 00406 * 00407 IF( M.LT.0 ) THEN 00408 INFO = -1 00409 ELSE IF( M.NE.N .AND. ISYM.NE.1 ) THEN 00410 INFO = -1 00411 ELSE IF( N.LT.0 ) THEN 00412 INFO = -2 00413 ELSE IF( IDIST.EQ.-1 ) THEN 00414 INFO = -3 00415 ELSE IF( ISYM.EQ.-1 ) THEN 00416 INFO = -5 00417 ELSE IF( ABS( MODE ).GT.6 ) THEN 00418 INFO = -7 00419 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE ) 00420 $ THEN 00421 INFO = -8 00422 ELSE IF( KL.LT.0 ) THEN 00423 INFO = -10 00424 ELSE IF( KU.LT.0 .OR. ( ISYM.NE.1 .AND. KL.NE.KU ) ) THEN 00425 INFO = -11 00426 ELSE IF( IPACK.EQ.-1 .OR. ( ISYMPK.EQ.1 .AND. ISYM.EQ.1 ) .OR. 00427 $ ( ISYMPK.EQ.2 .AND. ISYM.EQ.1 .AND. KL.GT.0 ) .OR. 00428 $ ( ISYMPK.EQ.3 .AND. ISYM.EQ.1 .AND. KU.GT.0 ) .OR. 00429 $ ( ISYMPK.NE.0 .AND. M.NE.N ) ) THEN 00430 INFO = -12 00431 ELSE IF( LDA.LT.MAX( 1, MINLDA ) ) THEN 00432 INFO = -14 00433 END IF 00434 * 00435 IF( INFO.NE.0 ) THEN 00436 CALL XERBLA( 'CLATMS', -INFO ) 00437 RETURN 00438 END IF 00439 * 00440 * Initialize random number generator 00441 * 00442 DO 10 I = 1, 4 00443 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 ) 00444 10 CONTINUE 00445 * 00446 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 00447 $ ISEED( 4 ) = ISEED( 4 ) + 1 00448 * 00449 * 2) Set up D if indicated. 00450 * 00451 * Compute D according to COND and MODE 00452 * 00453 CALL SLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, MNMIN, IINFO ) 00454 IF( IINFO.NE.0 ) THEN 00455 INFO = 1 00456 RETURN 00457 END IF 00458 * 00459 * Choose Top-Down if D is (apparently) increasing, 00460 * Bottom-Up if D is (apparently) decreasing. 00461 * 00462 IF( ABS( D( 1 ) ).LE.ABS( D( MNMIN ) ) ) THEN 00463 TOPDWN = .TRUE. 00464 ELSE 00465 TOPDWN = .FALSE. 00466 END IF 00467 * 00468 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN 00469 * 00470 * Scale by DMAX 00471 * 00472 TEMP = ABS( D( 1 ) ) 00473 DO 20 I = 2, MNMIN 00474 TEMP = MAX( TEMP, ABS( D( I ) ) ) 00475 20 CONTINUE 00476 * 00477 IF( TEMP.GT.ZERO ) THEN 00478 ALPHA = DMAX / TEMP 00479 ELSE 00480 INFO = 2 00481 RETURN 00482 END IF 00483 * 00484 CALL SSCAL( MNMIN, ALPHA, D, 1 ) 00485 * 00486 END IF 00487 * 00488 CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA ) 00489 * 00490 * 3) Generate Banded Matrix using Givens rotations. 00491 * Also the special case of UUB=LLB=0 00492 * 00493 * Compute Addressing constants to cover all 00494 * storage formats. Whether GE, HE, SY, GB, HB, or SB, 00495 * upper or lower triangle or both, 00496 * the (i,j)-th element is in 00497 * A( i - ISKEW*j + IOFFST, j ) 00498 * 00499 IF( IPACK.GT.4 ) THEN 00500 ILDA = LDA - 1 00501 ISKEW = 1 00502 IF( IPACK.GT.5 ) THEN 00503 IOFFST = UUB + 1 00504 ELSE 00505 IOFFST = 1 00506 END IF 00507 ELSE 00508 ILDA = LDA 00509 ISKEW = 0 00510 IOFFST = 0 00511 END IF 00512 * 00513 * IPACKG is the format that the matrix is generated in. If this is 00514 * different from IPACK, then the matrix must be repacked at the 00515 * end. It also signals how to compute the norm, for scaling. 00516 * 00517 IPACKG = 0 00518 * 00519 * Diagonal Matrix -- We are done, unless it 00520 * is to be stored HP/SP/PP/TP (PACK='R' or 'C') 00521 * 00522 IF( LLB.EQ.0 .AND. UUB.EQ.0 ) THEN 00523 DO 30 J = 1, MNMIN 00524 A( ( 1-ISKEW )*J+IOFFST, J ) = CMPLX( D( J ) ) 00525 30 CONTINUE 00526 * 00527 IF( IPACK.LE.2 .OR. IPACK.GE.5 ) 00528 $ IPACKG = IPACK 00529 * 00530 ELSE IF( GIVENS ) THEN 00531 * 00532 * Check whether to use Givens rotations, 00533 * Householder transformations, or nothing. 00534 * 00535 IF( ISYM.EQ.1 ) THEN 00536 * 00537 * Non-symmetric -- A = U D V 00538 * 00539 IF( IPACK.GT.4 ) THEN 00540 IPACKG = IPACK 00541 ELSE 00542 IPACKG = 0 00543 END IF 00544 * 00545 DO 40 J = 1, MNMIN 00546 A( ( 1-ISKEW )*J+IOFFST, J ) = CMPLX( D( J ) ) 00547 40 CONTINUE 00548 * 00549 IF( TOPDWN ) THEN 00550 JKL = 0 00551 DO 70 JKU = 1, UUB 00552 * 00553 * Transform from bandwidth JKL, JKU-1 to JKL, JKU 00554 * 00555 * Last row actually rotated is M 00556 * Last column actually rotated is MIN( M+JKU, N ) 00557 * 00558 DO 60 JR = 1, MIN( M+JKU, N ) + JKL - 1 00559 EXTRA = CZERO 00560 ANGLE = TWOPI*SLARND( 1, ISEED ) 00561 C = COS( ANGLE )*CLARND( 5, ISEED ) 00562 S = SIN( ANGLE )*CLARND( 5, ISEED ) 00563 ICOL = MAX( 1, JR-JKL ) 00564 IF( JR.LT.M ) THEN 00565 IL = MIN( N, JR+JKU ) + 1 - ICOL 00566 CALL CLAROT( .TRUE., JR.GT.JKL, .FALSE., IL, C, 00567 $ S, A( JR-ISKEW*ICOL+IOFFST, ICOL ), 00568 $ ILDA, EXTRA, DUMMY ) 00569 END IF 00570 * 00571 * Chase "EXTRA" back up 00572 * 00573 IR = JR 00574 IC = ICOL 00575 DO 50 JCH = JR - JKL, 1, -JKL - JKU 00576 IF( IR.LT.M ) THEN 00577 CALL CLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST, 00578 $ IC+1 ), EXTRA, REALC, S, DUMMY ) 00579 DUMMY = CLARND( 5, ISEED ) 00580 C = CONJG( REALC*DUMMY ) 00581 S = CONJG( -S*DUMMY ) 00582 END IF 00583 IROW = MAX( 1, JCH-JKU ) 00584 IL = IR + 2 - IROW 00585 CTEMP = CZERO 00586 ILTEMP = JCH.GT.JKU 00587 CALL CLAROT( .FALSE., ILTEMP, .TRUE., IL, C, S, 00588 $ A( IROW-ISKEW*IC+IOFFST, IC ), 00589 $ ILDA, CTEMP, EXTRA ) 00590 IF( ILTEMP ) THEN 00591 CALL CLARTG( A( IROW+1-ISKEW*( IC+1 )+IOFFST, 00592 $ IC+1 ), CTEMP, REALC, S, DUMMY ) 00593 DUMMY = CLARND( 5, ISEED ) 00594 C = CONJG( REALC*DUMMY ) 00595 S = CONJG( -S*DUMMY ) 00596 * 00597 ICOL = MAX( 1, JCH-JKU-JKL ) 00598 IL = IC + 2 - ICOL 00599 EXTRA = CZERO 00600 CALL CLAROT( .TRUE., JCH.GT.JKU+JKL, .TRUE., 00601 $ IL, C, S, A( IROW-ISKEW*ICOL+ 00602 $ IOFFST, ICOL ), ILDA, EXTRA, 00603 $ CTEMP ) 00604 IC = ICOL 00605 IR = IROW 00606 END IF 00607 50 CONTINUE 00608 60 CONTINUE 00609 70 CONTINUE 00610 * 00611 JKU = UUB 00612 DO 100 JKL = 1, LLB 00613 * 00614 * Transform from bandwidth JKL-1, JKU to JKL, JKU 00615 * 00616 DO 90 JC = 1, MIN( N+JKL, M ) + JKU - 1 00617 EXTRA = CZERO 00618 ANGLE = TWOPI*SLARND( 1, ISEED ) 00619 C = COS( ANGLE )*CLARND( 5, ISEED ) 00620 S = SIN( ANGLE )*CLARND( 5, ISEED ) 00621 IROW = MAX( 1, JC-JKU ) 00622 IF( JC.LT.N ) THEN 00623 IL = MIN( M, JC+JKL ) + 1 - IROW 00624 CALL CLAROT( .FALSE., JC.GT.JKU, .FALSE., IL, C, 00625 $ S, A( IROW-ISKEW*JC+IOFFST, JC ), 00626 $ ILDA, EXTRA, DUMMY ) 00627 END IF 00628 * 00629 * Chase "EXTRA" back up 00630 * 00631 IC = JC 00632 IR = IROW 00633 DO 80 JCH = JC - JKU, 1, -JKL - JKU 00634 IF( IC.LT.N ) THEN 00635 CALL CLARTG( A( IR+1-ISKEW*( IC+1 )+IOFFST, 00636 $ IC+1 ), EXTRA, REALC, S, DUMMY ) 00637 DUMMY = CLARND( 5, ISEED ) 00638 C = CONJG( REALC*DUMMY ) 00639 S = CONJG( -S*DUMMY ) 00640 END IF 00641 ICOL = MAX( 1, JCH-JKL ) 00642 IL = IC + 2 - ICOL 00643 CTEMP = CZERO 00644 ILTEMP = JCH.GT.JKL 00645 CALL CLAROT( .TRUE., ILTEMP, .TRUE., IL, C, S, 00646 $ A( IR-ISKEW*ICOL+IOFFST, ICOL ), 00647 $ ILDA, CTEMP, EXTRA ) 00648 IF( ILTEMP ) THEN 00649 CALL CLARTG( A( IR+1-ISKEW*( ICOL+1 )+IOFFST, 00650 $ ICOL+1 ), CTEMP, REALC, S, 00651 $ DUMMY ) 00652 DUMMY = CLARND( 5, ISEED ) 00653 C = CONJG( REALC*DUMMY ) 00654 S = CONJG( -S*DUMMY ) 00655 IROW = MAX( 1, JCH-JKL-JKU ) 00656 IL = IR + 2 - IROW 00657 EXTRA = CZERO 00658 CALL CLAROT( .FALSE., JCH.GT.JKL+JKU, .TRUE., 00659 $ IL, C, S, A( IROW-ISKEW*ICOL+ 00660 $ IOFFST, ICOL ), ILDA, EXTRA, 00661 $ CTEMP ) 00662 IC = ICOL 00663 IR = IROW 00664 END IF 00665 80 CONTINUE 00666 90 CONTINUE 00667 100 CONTINUE 00668 * 00669 ELSE 00670 * 00671 * Bottom-Up -- Start at the bottom right. 00672 * 00673 JKL = 0 00674 DO 130 JKU = 1, UUB 00675 * 00676 * Transform from bandwidth JKL, JKU-1 to JKL, JKU 00677 * 00678 * First row actually rotated is M 00679 * First column actually rotated is MIN( M+JKU, N ) 00680 * 00681 IENDCH = MIN( M, N+JKL ) - 1 00682 DO 120 JC = MIN( M+JKU, N ) - 1, 1 - JKL, -1 00683 EXTRA = CZERO 00684 ANGLE = TWOPI*SLARND( 1, ISEED ) 00685 C = COS( ANGLE )*CLARND( 5, ISEED ) 00686 S = SIN( ANGLE )*CLARND( 5, ISEED ) 00687 IROW = MAX( 1, JC-JKU+1 ) 00688 IF( JC.GT.0 ) THEN 00689 IL = MIN( M, JC+JKL+1 ) + 1 - IROW 00690 CALL CLAROT( .FALSE., .FALSE., JC+JKL.LT.M, IL, 00691 $ C, S, A( IROW-ISKEW*JC+IOFFST, 00692 $ JC ), ILDA, DUMMY, EXTRA ) 00693 END IF 00694 * 00695 * Chase "EXTRA" back down 00696 * 00697 IC = JC 00698 DO 110 JCH = JC + JKL, IENDCH, JKL + JKU 00699 ILEXTR = IC.GT.0 00700 IF( ILEXTR ) THEN 00701 CALL CLARTG( A( JCH-ISKEW*IC+IOFFST, IC ), 00702 $ EXTRA, REALC, S, DUMMY ) 00703 DUMMY = CLARND( 5, ISEED ) 00704 C = REALC*DUMMY 00705 S = S*DUMMY 00706 END IF 00707 IC = MAX( 1, IC ) 00708 ICOL = MIN( N-1, JCH+JKU ) 00709 ILTEMP = JCH + JKU.LT.N 00710 CTEMP = CZERO 00711 CALL CLAROT( .TRUE., ILEXTR, ILTEMP, ICOL+2-IC, 00712 $ C, S, A( JCH-ISKEW*IC+IOFFST, IC ), 00713 $ ILDA, EXTRA, CTEMP ) 00714 IF( ILTEMP ) THEN 00715 CALL CLARTG( A( JCH-ISKEW*ICOL+IOFFST, 00716 $ ICOL ), CTEMP, REALC, S, DUMMY ) 00717 DUMMY = CLARND( 5, ISEED ) 00718 C = REALC*DUMMY 00719 S = S*DUMMY 00720 IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH 00721 EXTRA = CZERO 00722 CALL CLAROT( .FALSE., .TRUE., 00723 $ JCH+JKL+JKU.LE.IENDCH, IL, C, S, 00724 $ A( JCH-ISKEW*ICOL+IOFFST, 00725 $ ICOL ), ILDA, CTEMP, EXTRA ) 00726 IC = ICOL 00727 END IF 00728 110 CONTINUE 00729 120 CONTINUE 00730 130 CONTINUE 00731 * 00732 JKU = UUB 00733 DO 160 JKL = 1, LLB 00734 * 00735 * Transform from bandwidth JKL-1, JKU to JKL, JKU 00736 * 00737 * First row actually rotated is MIN( N+JKL, M ) 00738 * First column actually rotated is N 00739 * 00740 IENDCH = MIN( N, M+JKU ) - 1 00741 DO 150 JR = MIN( N+JKL, M ) - 1, 1 - JKU, -1 00742 EXTRA = CZERO 00743 ANGLE = TWOPI*SLARND( 1, ISEED ) 00744 C = COS( ANGLE )*CLARND( 5, ISEED ) 00745 S = SIN( ANGLE )*CLARND( 5, ISEED ) 00746 ICOL = MAX( 1, JR-JKL+1 ) 00747 IF( JR.GT.0 ) THEN 00748 IL = MIN( N, JR+JKU+1 ) + 1 - ICOL 00749 CALL CLAROT( .TRUE., .FALSE., JR+JKU.LT.N, IL, 00750 $ C, S, A( JR-ISKEW*ICOL+IOFFST, 00751 $ ICOL ), ILDA, DUMMY, EXTRA ) 00752 END IF 00753 * 00754 * Chase "EXTRA" back down 00755 * 00756 IR = JR 00757 DO 140 JCH = JR + JKU, IENDCH, JKL + JKU 00758 ILEXTR = IR.GT.0 00759 IF( ILEXTR ) THEN 00760 CALL CLARTG( A( IR-ISKEW*JCH+IOFFST, JCH ), 00761 $ EXTRA, REALC, S, DUMMY ) 00762 DUMMY = CLARND( 5, ISEED ) 00763 C = REALC*DUMMY 00764 S = S*DUMMY 00765 END IF 00766 IR = MAX( 1, IR ) 00767 IROW = MIN( M-1, JCH+JKL ) 00768 ILTEMP = JCH + JKL.LT.M 00769 CTEMP = CZERO 00770 CALL CLAROT( .FALSE., ILEXTR, ILTEMP, IROW+2-IR, 00771 $ C, S, A( IR-ISKEW*JCH+IOFFST, 00772 $ JCH ), ILDA, EXTRA, CTEMP ) 00773 IF( ILTEMP ) THEN 00774 CALL CLARTG( A( IROW-ISKEW*JCH+IOFFST, JCH ), 00775 $ CTEMP, REALC, S, DUMMY ) 00776 DUMMY = CLARND( 5, ISEED ) 00777 C = REALC*DUMMY 00778 S = S*DUMMY 00779 IL = MIN( IENDCH, JCH+JKL+JKU ) + 2 - JCH 00780 EXTRA = CZERO 00781 CALL CLAROT( .TRUE., .TRUE., 00782 $ JCH+JKL+JKU.LE.IENDCH, IL, C, S, 00783 $ A( IROW-ISKEW*JCH+IOFFST, JCH ), 00784 $ ILDA, CTEMP, EXTRA ) 00785 IR = IROW 00786 END IF 00787 140 CONTINUE 00788 150 CONTINUE 00789 160 CONTINUE 00790 * 00791 END IF 00792 * 00793 ELSE 00794 * 00795 * Symmetric -- A = U D U' 00796 * Hermitian -- A = U D U* 00797 * 00798 IPACKG = IPACK 00799 IOFFG = IOFFST 00800 * 00801 IF( TOPDWN ) THEN 00802 * 00803 * Top-Down -- Generate Upper triangle only 00804 * 00805 IF( IPACK.GE.5 ) THEN 00806 IPACKG = 6 00807 IOFFG = UUB + 1 00808 ELSE 00809 IPACKG = 1 00810 END IF 00811 * 00812 DO 170 J = 1, MNMIN 00813 A( ( 1-ISKEW )*J+IOFFG, J ) = CMPLX( D( J ) ) 00814 170 CONTINUE 00815 * 00816 DO 200 K = 1, UUB 00817 DO 190 JC = 1, N - 1 00818 IROW = MAX( 1, JC-K ) 00819 IL = MIN( JC+1, K+2 ) 00820 EXTRA = CZERO 00821 CTEMP = A( JC-ISKEW*( JC+1 )+IOFFG, JC+1 ) 00822 ANGLE = TWOPI*SLARND( 1, ISEED ) 00823 C = COS( ANGLE )*CLARND( 5, ISEED ) 00824 S = SIN( ANGLE )*CLARND( 5, ISEED ) 00825 IF( CSYM ) THEN 00826 CT = C 00827 ST = S 00828 ELSE 00829 CTEMP = CONJG( CTEMP ) 00830 CT = CONJG( C ) 00831 ST = CONJG( S ) 00832 END IF 00833 CALL CLAROT( .FALSE., JC.GT.K, .TRUE., IL, C, S, 00834 $ A( IROW-ISKEW*JC+IOFFG, JC ), ILDA, 00835 $ EXTRA, CTEMP ) 00836 CALL CLAROT( .TRUE., .TRUE., .FALSE., 00837 $ MIN( K, N-JC )+1, CT, ST, 00838 $ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA, 00839 $ CTEMP, DUMMY ) 00840 * 00841 * Chase EXTRA back up the matrix 00842 * 00843 ICOL = JC 00844 DO 180 JCH = JC - K, 1, -K 00845 CALL CLARTG( A( JCH+1-ISKEW*( ICOL+1 )+IOFFG, 00846 $ ICOL+1 ), EXTRA, REALC, S, DUMMY ) 00847 DUMMY = CLARND( 5, ISEED ) 00848 C = CONJG( REALC*DUMMY ) 00849 S = CONJG( -S*DUMMY ) 00850 CTEMP = A( JCH-ISKEW*( JCH+1 )+IOFFG, JCH+1 ) 00851 IF( CSYM ) THEN 00852 CT = C 00853 ST = S 00854 ELSE 00855 CTEMP = CONJG( CTEMP ) 00856 CT = CONJG( C ) 00857 ST = CONJG( S ) 00858 END IF 00859 CALL CLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S, 00860 $ A( ( 1-ISKEW )*JCH+IOFFG, JCH ), 00861 $ ILDA, CTEMP, EXTRA ) 00862 IROW = MAX( 1, JCH-K ) 00863 IL = MIN( JCH+1, K+2 ) 00864 EXTRA = CZERO 00865 CALL CLAROT( .FALSE., JCH.GT.K, .TRUE., IL, CT, 00866 $ ST, A( IROW-ISKEW*JCH+IOFFG, JCH ), 00867 $ ILDA, EXTRA, CTEMP ) 00868 ICOL = JCH 00869 180 CONTINUE 00870 190 CONTINUE 00871 200 CONTINUE 00872 * 00873 * If we need lower triangle, copy from upper. Note that 00874 * the order of copying is chosen to work for 'q' -> 'b' 00875 * 00876 IF( IPACK.NE.IPACKG .AND. IPACK.NE.3 ) THEN 00877 DO 230 JC = 1, N 00878 IROW = IOFFST - ISKEW*JC 00879 IF( CSYM ) THEN 00880 DO 210 JR = JC, MIN( N, JC+UUB ) 00881 A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR ) 00882 210 CONTINUE 00883 ELSE 00884 DO 220 JR = JC, MIN( N, JC+UUB ) 00885 A( JR+IROW, JC ) = CONJG( A( JC-ISKEW*JR+ 00886 $ IOFFG, JR ) ) 00887 220 CONTINUE 00888 END IF 00889 230 CONTINUE 00890 IF( IPACK.EQ.5 ) THEN 00891 DO 250 JC = N - UUB + 1, N 00892 DO 240 JR = N + 2 - JC, UUB + 1 00893 A( JR, JC ) = CZERO 00894 240 CONTINUE 00895 250 CONTINUE 00896 END IF 00897 IF( IPACKG.EQ.6 ) THEN 00898 IPACKG = IPACK 00899 ELSE 00900 IPACKG = 0 00901 END IF 00902 END IF 00903 ELSE 00904 * 00905 * Bottom-Up -- Generate Lower triangle only 00906 * 00907 IF( IPACK.GE.5 ) THEN 00908 IPACKG = 5 00909 IF( IPACK.EQ.6 ) 00910 $ IOFFG = 1 00911 ELSE 00912 IPACKG = 2 00913 END IF 00914 * 00915 DO 260 J = 1, MNMIN 00916 A( ( 1-ISKEW )*J+IOFFG, J ) = CMPLX( D( J ) ) 00917 260 CONTINUE 00918 * 00919 DO 290 K = 1, UUB 00920 DO 280 JC = N - 1, 1, -1 00921 IL = MIN( N+1-JC, K+2 ) 00922 EXTRA = CZERO 00923 CTEMP = A( 1+( 1-ISKEW )*JC+IOFFG, JC ) 00924 ANGLE = TWOPI*SLARND( 1, ISEED ) 00925 C = COS( ANGLE )*CLARND( 5, ISEED ) 00926 S = SIN( ANGLE )*CLARND( 5, ISEED ) 00927 IF( CSYM ) THEN 00928 CT = C 00929 ST = S 00930 ELSE 00931 CTEMP = CONJG( CTEMP ) 00932 CT = CONJG( C ) 00933 ST = CONJG( S ) 00934 END IF 00935 CALL CLAROT( .FALSE., .TRUE., N-JC.GT.K, IL, C, S, 00936 $ A( ( 1-ISKEW )*JC+IOFFG, JC ), ILDA, 00937 $ CTEMP, EXTRA ) 00938 ICOL = MAX( 1, JC-K+1 ) 00939 CALL CLAROT( .TRUE., .FALSE., .TRUE., JC+2-ICOL, 00940 $ CT, ST, A( JC-ISKEW*ICOL+IOFFG, 00941 $ ICOL ), ILDA, DUMMY, CTEMP ) 00942 * 00943 * Chase EXTRA back down the matrix 00944 * 00945 ICOL = JC 00946 DO 270 JCH = JC + K, N - 1, K 00947 CALL CLARTG( A( JCH-ISKEW*ICOL+IOFFG, ICOL ), 00948 $ EXTRA, REALC, S, DUMMY ) 00949 DUMMY = CLARND( 5, ISEED ) 00950 C = REALC*DUMMY 00951 S = S*DUMMY 00952 CTEMP = A( 1+( 1-ISKEW )*JCH+IOFFG, JCH ) 00953 IF( CSYM ) THEN 00954 CT = C 00955 ST = S 00956 ELSE 00957 CTEMP = CONJG( CTEMP ) 00958 CT = CONJG( C ) 00959 ST = CONJG( S ) 00960 END IF 00961 CALL CLAROT( .TRUE., .TRUE., .TRUE., K+2, C, S, 00962 $ A( JCH-ISKEW*ICOL+IOFFG, ICOL ), 00963 $ ILDA, EXTRA, CTEMP ) 00964 IL = MIN( N+1-JCH, K+2 ) 00965 EXTRA = CZERO 00966 CALL CLAROT( .FALSE., .TRUE., N-JCH.GT.K, IL, 00967 $ CT, ST, A( ( 1-ISKEW )*JCH+IOFFG, 00968 $ JCH ), ILDA, CTEMP, EXTRA ) 00969 ICOL = JCH 00970 270 CONTINUE 00971 280 CONTINUE 00972 290 CONTINUE 00973 * 00974 * If we need upper triangle, copy from lower. Note that 00975 * the order of copying is chosen to work for 'b' -> 'q' 00976 * 00977 IF( IPACK.NE.IPACKG .AND. IPACK.NE.4 ) THEN 00978 DO 320 JC = N, 1, -1 00979 IROW = IOFFST - ISKEW*JC 00980 IF( CSYM ) THEN 00981 DO 300 JR = JC, MAX( 1, JC-UUB ), -1 00982 A( JR+IROW, JC ) = A( JC-ISKEW*JR+IOFFG, JR ) 00983 300 CONTINUE 00984 ELSE 00985 DO 310 JR = JC, MAX( 1, JC-UUB ), -1 00986 A( JR+IROW, JC ) = CONJG( A( JC-ISKEW*JR+ 00987 $ IOFFG, JR ) ) 00988 310 CONTINUE 00989 END IF 00990 320 CONTINUE 00991 IF( IPACK.EQ.6 ) THEN 00992 DO 340 JC = 1, UUB 00993 DO 330 JR = 1, UUB + 1 - JC 00994 A( JR, JC ) = CZERO 00995 330 CONTINUE 00996 340 CONTINUE 00997 END IF 00998 IF( IPACKG.EQ.5 ) THEN 00999 IPACKG = IPACK 01000 ELSE 01001 IPACKG = 0 01002 END IF 01003 END IF 01004 END IF 01005 * 01006 * Ensure that the diagonal is real if Hermitian 01007 * 01008 IF( .NOT.CSYM ) THEN 01009 DO 350 JC = 1, N 01010 IROW = IOFFST + ( 1-ISKEW )*JC 01011 A( IROW, JC ) = CMPLX( REAL( A( IROW, JC ) ) ) 01012 350 CONTINUE 01013 END IF 01014 * 01015 END IF 01016 * 01017 ELSE 01018 * 01019 * 4) Generate Banded Matrix by first 01020 * Rotating by random Unitary matrices, 01021 * then reducing the bandwidth using Householder 01022 * transformations. 01023 * 01024 * Note: we should get here only if LDA .ge. N 01025 * 01026 IF( ISYM.EQ.1 ) THEN 01027 * 01028 * Non-symmetric -- A = U D V 01029 * 01030 CALL CLAGGE( MR, NC, LLB, UUB, D, A, LDA, ISEED, WORK, 01031 $ IINFO ) 01032 ELSE 01033 * 01034 * Symmetric -- A = U D U' or 01035 * Hermitian -- A = U D U* 01036 * 01037 IF( CSYM ) THEN 01038 CALL CLAGSY( M, LLB, D, A, LDA, ISEED, WORK, IINFO ) 01039 ELSE 01040 CALL CLAGHE( M, LLB, D, A, LDA, ISEED, WORK, IINFO ) 01041 END IF 01042 END IF 01043 * 01044 IF( IINFO.NE.0 ) THEN 01045 INFO = 3 01046 RETURN 01047 END IF 01048 END IF 01049 * 01050 * 5) Pack the matrix 01051 * 01052 IF( IPACK.NE.IPACKG ) THEN 01053 IF( IPACK.EQ.1 ) THEN 01054 * 01055 * 'U' -- Upper triangular, not packed 01056 * 01057 DO 370 J = 1, M 01058 DO 360 I = J + 1, M 01059 A( I, J ) = CZERO 01060 360 CONTINUE 01061 370 CONTINUE 01062 * 01063 ELSE IF( IPACK.EQ.2 ) THEN 01064 * 01065 * 'L' -- Lower triangular, not packed 01066 * 01067 DO 390 J = 2, M 01068 DO 380 I = 1, J - 1 01069 A( I, J ) = CZERO 01070 380 CONTINUE 01071 390 CONTINUE 01072 * 01073 ELSE IF( IPACK.EQ.3 ) THEN 01074 * 01075 * 'C' -- Upper triangle packed Columnwise. 01076 * 01077 ICOL = 1 01078 IROW = 0 01079 DO 410 J = 1, M 01080 DO 400 I = 1, J 01081 IROW = IROW + 1 01082 IF( IROW.GT.LDA ) THEN 01083 IROW = 1 01084 ICOL = ICOL + 1 01085 END IF 01086 A( IROW, ICOL ) = A( I, J ) 01087 400 CONTINUE 01088 410 CONTINUE 01089 * 01090 ELSE IF( IPACK.EQ.4 ) THEN 01091 * 01092 * 'R' -- Lower triangle packed Columnwise. 01093 * 01094 ICOL = 1 01095 IROW = 0 01096 DO 430 J = 1, M 01097 DO 420 I = J, M 01098 IROW = IROW + 1 01099 IF( IROW.GT.LDA ) THEN 01100 IROW = 1 01101 ICOL = ICOL + 1 01102 END IF 01103 A( IROW, ICOL ) = A( I, J ) 01104 420 CONTINUE 01105 430 CONTINUE 01106 * 01107 ELSE IF( IPACK.GE.5 ) THEN 01108 * 01109 * 'B' -- The lower triangle is packed as a band matrix. 01110 * 'Q' -- The upper triangle is packed as a band matrix. 01111 * 'Z' -- The whole matrix is packed as a band matrix. 01112 * 01113 IF( IPACK.EQ.5 ) 01114 $ UUB = 0 01115 IF( IPACK.EQ.6 ) 01116 $ LLB = 0 01117 * 01118 DO 450 J = 1, UUB 01119 DO 440 I = MIN( J+LLB, M ), 1, -1 01120 A( I-J+UUB+1, J ) = A( I, J ) 01121 440 CONTINUE 01122 450 CONTINUE 01123 * 01124 DO 470 J = UUB + 2, N 01125 DO 460 I = J - UUB, MIN( J+LLB, M ) 01126 A( I-J+UUB+1, J ) = A( I, J ) 01127 460 CONTINUE 01128 470 CONTINUE 01129 END IF 01130 * 01131 * If packed, zero out extraneous elements. 01132 * 01133 * Symmetric/Triangular Packed -- 01134 * zero out everything after A(IROW,ICOL) 01135 * 01136 IF( IPACK.EQ.3 .OR. IPACK.EQ.4 ) THEN 01137 DO 490 JC = ICOL, M 01138 DO 480 JR = IROW + 1, LDA 01139 A( JR, JC ) = CZERO 01140 480 CONTINUE 01141 IROW = 0 01142 490 CONTINUE 01143 * 01144 ELSE IF( IPACK.GE.5 ) THEN 01145 * 01146 * Packed Band -- 01147 * 1st row is now in A( UUB+2-j, j), zero above it 01148 * m-th row is now in A( M+UUB-j,j), zero below it 01149 * last non-zero diagonal is now in A( UUB+LLB+1,j ), 01150 * zero below it, too. 01151 * 01152 IR1 = UUB + LLB + 2 01153 IR2 = UUB + M + 2 01154 DO 520 JC = 1, N 01155 DO 500 JR = 1, UUB + 1 - JC 01156 A( JR, JC ) = CZERO 01157 500 CONTINUE 01158 DO 510 JR = MAX( 1, MIN( IR1, IR2-JC ) ), LDA 01159 A( JR, JC ) = CZERO 01160 510 CONTINUE 01161 520 CONTINUE 01162 END IF 01163 END IF 01164 * 01165 RETURN 01166 * 01167 * End of CLATMS 01168 * 01169 END
