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