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