![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLATME 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 DLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 00012 * RSIGN, 00013 * UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 00014 * A, 00015 * LDA, WORK, INFO ) 00016 * 00017 * .. Scalar Arguments .. 00018 * CHARACTER DIST, RSIGN, SIM, UPPER 00019 * INTEGER INFO, KL, KU, LDA, MODE, MODES, N 00020 * DOUBLE PRECISION ANORM, COND, CONDS, DMAX 00021 * .. 00022 * .. Array Arguments .. 00023 * CHARACTER EI( * ) 00024 * INTEGER ISEED( 4 ) 00025 * DOUBLE PRECISION A( LDA, * ), D( * ), DS( * ), WORK( * ) 00026 * .. 00027 * 00028 * 00029 *> \par Purpose: 00030 * ============= 00031 *> 00032 *> \verbatim 00033 *> 00034 *> DLATME generates random non-symmetric square matrices with 00035 *> specified eigenvalues for testing LAPACK programs. 00036 *> 00037 *> DLATME operates by applying the following sequence of 00038 *> operations: 00039 *> 00040 *> 1. Set the diagonal to D, where D may be input or 00041 *> computed according to MODE, COND, DMAX, and RSIGN 00042 *> as described below. 00043 *> 00044 *> 2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R', 00045 *> or MODE=5), certain pairs of adjacent elements of D are 00046 *> interpreted as the real and complex parts of a complex 00047 *> conjugate pair; A thus becomes block diagonal, with 1x1 00048 *> and 2x2 blocks. 00049 *> 00050 *> 3. If UPPER='T', the upper triangle of A is set to random values 00051 *> out of distribution DIST. 00052 *> 00053 *> 4. If SIM='T', A is multiplied on the left by a random matrix 00054 *> X, whose singular values are specified by DS, MODES, and 00055 *> CONDS, and on the right by X inverse. 00056 *> 00057 *> 5. If KL < N-1, the lower bandwidth is reduced to KL using 00058 *> Householder transformations. If KU < N-1, the upper 00059 *> bandwidth is reduced to KU. 00060 *> 00061 *> 6. If ANORM is not negative, the matrix is scaled to have 00062 *> maximum-element-norm ANORM. 00063 *> 00064 *> (Note: since the matrix cannot be reduced beyond Hessenberg form, 00065 *> no packing options are available.) 00066 *> \endverbatim 00067 * 00068 * Arguments: 00069 * ========== 00070 * 00071 *> \param[in] N 00072 *> \verbatim 00073 *> N is INTEGER 00074 *> The number of columns (or rows) of A. Not modified. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] DIST 00078 *> \verbatim 00079 *> DIST is CHARACTER*1 00080 *> On entry, DIST specifies the type of distribution to be used 00081 *> to generate the random eigen-/singular values, and for the 00082 *> upper triangle (see UPPER). 00083 *> 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform ) 00084 *> 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric ) 00085 *> 'N' => NORMAL( 0, 1 ) ( 'N' for normal ) 00086 *> Not modified. 00087 *> \endverbatim 00088 *> 00089 *> \param[in,out] ISEED 00090 *> \verbatim 00091 *> ISEED is INTEGER array, dimension ( 4 ) 00092 *> On entry ISEED specifies the seed of the random number 00093 *> generator. They should lie between 0 and 4095 inclusive, 00094 *> and ISEED(4) should be odd. The random number generator 00095 *> uses a linear congruential sequence limited to small 00096 *> integers, and so should produce machine independent 00097 *> random numbers. The values of ISEED are changed on 00098 *> exit, and can be used in the next call to DLATME 00099 *> to continue the same random number sequence. 00100 *> Changed on exit. 00101 *> \endverbatim 00102 *> 00103 *> \param[in,out] D 00104 *> \verbatim 00105 *> D is DOUBLE PRECISION array, dimension ( N ) 00106 *> This array is used to specify the eigenvalues of A. If 00107 *> MODE=0, then D is assumed to contain the eigenvalues (but 00108 *> see the description of EI), otherwise they will be 00109 *> computed according to MODE, COND, DMAX, and RSIGN and 00110 *> placed in D. 00111 *> Modified if MODE is nonzero. 00112 *> \endverbatim 00113 *> 00114 *> \param[in] MODE 00115 *> \verbatim 00116 *> MODE is INTEGER 00117 *> On entry this describes how the eigenvalues are to 00118 *> be specified: 00119 *> MODE = 0 means use D (with EI) as input 00120 *> MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND 00121 *> MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND 00122 *> MODE = 3 sets D(I)=COND**(-(I-1)/(N-1)) 00123 *> MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND) 00124 *> MODE = 5 sets D to random numbers in the range 00125 *> ( 1/COND , 1 ) such that their logarithms 00126 *> are uniformly distributed. Each odd-even pair 00127 *> of elements will be either used as two real 00128 *> eigenvalues or as the real and imaginary part 00129 *> of a complex conjugate pair of eigenvalues; 00130 *> the choice of which is done is random, with 00131 *> 50-50 probability, for each pair. 00132 *> MODE = 6 set D to random numbers from same distribution 00133 *> as the rest of the matrix. 00134 *> MODE < 0 has the same meaning as ABS(MODE), except that 00135 *> the order of the elements of D is reversed. 00136 *> Thus if MODE is between 1 and 4, D has entries ranging 00137 *> from 1 to 1/COND, if between -1 and -4, D has entries 00138 *> ranging from 1/COND to 1, 00139 *> Not modified. 00140 *> \endverbatim 00141 *> 00142 *> \param[in] COND 00143 *> \verbatim 00144 *> COND is DOUBLE PRECISION 00145 *> On entry, this is used as described under MODE above. 00146 *> If used, it must be >= 1. Not modified. 00147 *> \endverbatim 00148 *> 00149 *> \param[in] DMAX 00150 *> \verbatim 00151 *> DMAX is DOUBLE PRECISION 00152 *> If MODE is neither -6, 0 nor 6, the contents of D, as 00153 *> computed according to MODE and COND, will be scaled by 00154 *> DMAX / max(abs(D(i))). Note that DMAX need not be 00155 *> positive: if DMAX is negative (or zero), D will be 00156 *> scaled by a negative number (or zero). 00157 *> Not modified. 00158 *> \endverbatim 00159 *> 00160 *> \param[in] EI 00161 *> \verbatim 00162 *> EI is CHARACTER*1 array, dimension ( N ) 00163 *> If MODE is 0, and EI(1) is not ' ' (space character), 00164 *> this array specifies which elements of D (on input) are 00165 *> real eigenvalues and which are the real and imaginary parts 00166 *> of a complex conjugate pair of eigenvalues. The elements 00167 *> of EI may then only have the values 'R' and 'I'. If 00168 *> EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is 00169 *> CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex 00170 *> conjugate thereof. If EI(j)=EI(j+1)='R', then the j-th 00171 *> eigenvalue is D(j) (i.e., real). EI(1) may not be 'I', 00172 *> nor may two adjacent elements of EI both have the value 'I'. 00173 *> If MODE is not 0, then EI is ignored. If MODE is 0 and 00174 *> EI(1)=' ', then the eigenvalues will all be real. 00175 *> Not modified. 00176 *> \endverbatim 00177 *> 00178 *> \param[in] RSIGN 00179 *> \verbatim 00180 *> RSIGN is CHARACTER*1 00181 *> If MODE is not 0, 6, or -6, and RSIGN='T', then the 00182 *> elements of D, as computed according to MODE and COND, will 00183 *> be multiplied by a random sign (+1 or -1). If RSIGN='F', 00184 *> they will not be. RSIGN may only have the values 'T' or 00185 *> 'F'. 00186 *> Not modified. 00187 *> \endverbatim 00188 *> 00189 *> \param[in] UPPER 00190 *> \verbatim 00191 *> UPPER is CHARACTER*1 00192 *> If UPPER='T', then the elements of A above the diagonal 00193 *> (and above the 2x2 diagonal blocks, if A has complex 00194 *> eigenvalues) will be set to random numbers out of DIST. 00195 *> If UPPER='F', they will not. UPPER may only have the 00196 *> values 'T' or 'F'. 00197 *> Not modified. 00198 *> \endverbatim 00199 *> 00200 *> \param[in] SIM 00201 *> \verbatim 00202 *> SIM is CHARACTER*1 00203 *> If SIM='T', then A will be operated on by a "similarity 00204 *> transform", i.e., multiplied on the left by a matrix X and 00205 *> on the right by X inverse. X = U S V, where U and V are 00206 *> random unitary matrices and S is a (diagonal) matrix of 00207 *> singular values specified by DS, MODES, and CONDS. If 00208 *> SIM='F', then A will not be transformed. 00209 *> Not modified. 00210 *> \endverbatim 00211 *> 00212 *> \param[in,out] DS 00213 *> \verbatim 00214 *> DS is DOUBLE PRECISION array, dimension ( N ) 00215 *> This array is used to specify the singular values of X, 00216 *> in the same way that D specifies the eigenvalues of A. 00217 *> If MODE=0, the DS contains the singular values, which 00218 *> may not be zero. 00219 *> Modified if MODE is nonzero. 00220 *> \endverbatim 00221 *> 00222 *> \param[in] MODES 00223 *> \verbatim 00224 *> MODES is INTEGER 00225 *> \endverbatim 00226 *> 00227 *> \param[in] CONDS 00228 *> \verbatim 00229 *> CONDS is DOUBLE PRECISION 00230 *> Same as MODE and COND, but for specifying the diagonal 00231 *> of S. MODES=-6 and +6 are not allowed (since they would 00232 *> result in randomly ill-conditioned eigenvalues.) 00233 *> \endverbatim 00234 *> 00235 *> \param[in] KL 00236 *> \verbatim 00237 *> KL is INTEGER 00238 *> This specifies the lower bandwidth of the matrix. KL=1 00239 *> specifies upper Hessenberg form. If KL is at least N-1, 00240 *> then A will have full lower bandwidth. KL must be at 00241 *> least 1. 00242 *> Not modified. 00243 *> \endverbatim 00244 *> 00245 *> \param[in] KU 00246 *> \verbatim 00247 *> KU is INTEGER 00248 *> This specifies the upper bandwidth of the matrix. KU=1 00249 *> specifies lower Hessenberg form. If KU is at least N-1, 00250 *> then A will have full upper bandwidth; if KU and KL 00251 *> are both at least N-1, then A will be dense. Only one of 00252 *> KU and KL may be less than N-1. KU must be at least 1. 00253 *> Not modified. 00254 *> \endverbatim 00255 *> 00256 *> \param[in] ANORM 00257 *> \verbatim 00258 *> ANORM is DOUBLE PRECISION 00259 *> If ANORM is not negative, then A will be scaled by a non- 00260 *> negative real number to make the maximum-element-norm of A 00261 *> to be ANORM. 00262 *> Not modified. 00263 *> \endverbatim 00264 *> 00265 *> \param[out] A 00266 *> \verbatim 00267 *> A is DOUBLE PRECISION array, dimension ( LDA, N ) 00268 *> On exit A is the desired test matrix. 00269 *> Modified. 00270 *> \endverbatim 00271 *> 00272 *> \param[in] LDA 00273 *> \verbatim 00274 *> LDA is INTEGER 00275 *> LDA specifies the first dimension of A as declared in the 00276 *> calling program. LDA must be at least N. 00277 *> Not modified. 00278 *> \endverbatim 00279 *> 00280 *> \param[out] WORK 00281 *> \verbatim 00282 *> WORK is DOUBLE PRECISION array, dimension ( 3*N ) 00283 *> Workspace. 00284 *> Modified. 00285 *> \endverbatim 00286 *> 00287 *> \param[out] INFO 00288 *> \verbatim 00289 *> INFO is INTEGER 00290 *> Error code. On exit, INFO will be set to one of the 00291 *> following values: 00292 *> 0 => normal return 00293 *> -1 => N negative 00294 *> -2 => DIST illegal string 00295 *> -5 => MODE not in range -6 to 6 00296 *> -6 => COND less than 1.0, and MODE neither -6, 0 nor 6 00297 *> -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or 00298 *> two adjacent elements of EI are 'I'. 00299 *> -9 => RSIGN is not 'T' or 'F' 00300 *> -10 => UPPER is not 'T' or 'F' 00301 *> -11 => SIM is not 'T' or 'F' 00302 *> -12 => MODES=0 and DS has a zero singular value. 00303 *> -13 => MODES is not in the range -5 to 5. 00304 *> -14 => MODES is nonzero and CONDS is less than 1. 00305 *> -15 => KL is less than 1. 00306 *> -16 => KU is less than 1, or KL and KU are both less than 00307 *> N-1. 00308 *> -19 => LDA is less than N. 00309 *> 1 => Error return from DLATM1 (computing D) 00310 *> 2 => Cannot scale to DMAX (max. eigenvalue is 0) 00311 *> 3 => Error return from DLATM1 (computing DS) 00312 *> 4 => Error return from DLARGE 00313 *> 5 => Zero singular value from DLATM1. 00314 *> \endverbatim 00315 * 00316 * Authors: 00317 * ======== 00318 * 00319 *> \author Univ. of Tennessee 00320 *> \author Univ. of California Berkeley 00321 *> \author Univ. of Colorado Denver 00322 *> \author NAG Ltd. 00323 * 00324 *> \date November 2011 00325 * 00326 *> \ingroup double_matgen 00327 * 00328 * ===================================================================== 00329 SUBROUTINE DLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI, 00330 $ RSIGN, 00331 $ UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, 00332 $ A, 00333 $ LDA, WORK, INFO ) 00334 * 00335 * -- LAPACK computational routine (version 3.4.0) -- 00336 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00338 * November 2011 00339 * 00340 * .. Scalar Arguments .. 00341 CHARACTER DIST, RSIGN, SIM, UPPER 00342 INTEGER INFO, KL, KU, LDA, MODE, MODES, N 00343 DOUBLE PRECISION ANORM, COND, CONDS, DMAX 00344 * .. 00345 * .. Array Arguments .. 00346 CHARACTER EI( * ) 00347 INTEGER ISEED( 4 ) 00348 DOUBLE PRECISION A( LDA, * ), D( * ), DS( * ), WORK( * ) 00349 * .. 00350 * 00351 * ===================================================================== 00352 * 00353 * .. Parameters .. 00354 DOUBLE PRECISION ZERO 00355 PARAMETER ( ZERO = 0.0D0 ) 00356 DOUBLE PRECISION ONE 00357 PARAMETER ( ONE = 1.0D0 ) 00358 DOUBLE PRECISION HALF 00359 PARAMETER ( HALF = 1.0D0 / 2.0D0 ) 00360 * .. 00361 * .. Local Scalars .. 00362 LOGICAL BADEI, BADS, USEEI 00363 INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN, 00364 $ ISIM, IUPPER, J, JC, JCR, JR 00365 DOUBLE PRECISION ALPHA, TAU, TEMP, XNORMS 00366 * .. 00367 * .. Local Arrays .. 00368 DOUBLE PRECISION TEMPA( 1 ) 00369 * .. 00370 * .. External Functions .. 00371 LOGICAL LSAME 00372 DOUBLE PRECISION DLANGE, DLARAN 00373 EXTERNAL LSAME, DLANGE, DLARAN 00374 * .. 00375 * .. External Subroutines .. 00376 EXTERNAL DCOPY, DGEMV, DGER, DLARFG, DLARGE, DLARNV, 00377 $ DLASET, DLATM1, DSCAL, XERBLA 00378 * .. 00379 * .. Intrinsic Functions .. 00380 INTRINSIC ABS, MAX, MOD 00381 * .. 00382 * .. Executable Statements .. 00383 * 00384 * 1) Decode and Test the input parameters. 00385 * Initialize flags & seed. 00386 * 00387 INFO = 0 00388 * 00389 * Quick return if possible 00390 * 00391 IF( N.EQ.0 ) 00392 $ RETURN 00393 * 00394 * Decode DIST 00395 * 00396 IF( LSAME( DIST, 'U' ) ) THEN 00397 IDIST = 1 00398 ELSE IF( LSAME( DIST, 'S' ) ) THEN 00399 IDIST = 2 00400 ELSE IF( LSAME( DIST, 'N' ) ) THEN 00401 IDIST = 3 00402 ELSE 00403 IDIST = -1 00404 END IF 00405 * 00406 * Check EI 00407 * 00408 USEEI = .TRUE. 00409 BADEI = .FALSE. 00410 IF( LSAME( EI( 1 ), ' ' ) .OR. MODE.NE.0 ) THEN 00411 USEEI = .FALSE. 00412 ELSE 00413 IF( LSAME( EI( 1 ), 'R' ) ) THEN 00414 DO 10 J = 2, N 00415 IF( LSAME( EI( J ), 'I' ) ) THEN 00416 IF( LSAME( EI( J-1 ), 'I' ) ) 00417 $ BADEI = .TRUE. 00418 ELSE 00419 IF( .NOT.LSAME( EI( J ), 'R' ) ) 00420 $ BADEI = .TRUE. 00421 END IF 00422 10 CONTINUE 00423 ELSE 00424 BADEI = .TRUE. 00425 END IF 00426 END IF 00427 * 00428 * Decode RSIGN 00429 * 00430 IF( LSAME( RSIGN, 'T' ) ) THEN 00431 IRSIGN = 1 00432 ELSE IF( LSAME( RSIGN, 'F' ) ) THEN 00433 IRSIGN = 0 00434 ELSE 00435 IRSIGN = -1 00436 END IF 00437 * 00438 * Decode UPPER 00439 * 00440 IF( LSAME( UPPER, 'T' ) ) THEN 00441 IUPPER = 1 00442 ELSE IF( LSAME( UPPER, 'F' ) ) THEN 00443 IUPPER = 0 00444 ELSE 00445 IUPPER = -1 00446 END IF 00447 * 00448 * Decode SIM 00449 * 00450 IF( LSAME( SIM, 'T' ) ) THEN 00451 ISIM = 1 00452 ELSE IF( LSAME( SIM, 'F' ) ) THEN 00453 ISIM = 0 00454 ELSE 00455 ISIM = -1 00456 END IF 00457 * 00458 * Check DS, if MODES=0 and ISIM=1 00459 * 00460 BADS = .FALSE. 00461 IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN 00462 DO 20 J = 1, N 00463 IF( DS( J ).EQ.ZERO ) 00464 $ BADS = .TRUE. 00465 20 CONTINUE 00466 END IF 00467 * 00468 * Set INFO if an error 00469 * 00470 IF( N.LT.0 ) THEN 00471 INFO = -1 00472 ELSE IF( IDIST.EQ.-1 ) THEN 00473 INFO = -2 00474 ELSE IF( ABS( MODE ).GT.6 ) THEN 00475 INFO = -5 00476 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE ) 00477 $ THEN 00478 INFO = -6 00479 ELSE IF( BADEI ) THEN 00480 INFO = -8 00481 ELSE IF( IRSIGN.EQ.-1 ) THEN 00482 INFO = -9 00483 ELSE IF( IUPPER.EQ.-1 ) THEN 00484 INFO = -10 00485 ELSE IF( ISIM.EQ.-1 ) THEN 00486 INFO = -11 00487 ELSE IF( BADS ) THEN 00488 INFO = -12 00489 ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN 00490 INFO = -13 00491 ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN 00492 INFO = -14 00493 ELSE IF( KL.LT.1 ) THEN 00494 INFO = -15 00495 ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN 00496 INFO = -16 00497 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00498 INFO = -19 00499 END IF 00500 * 00501 IF( INFO.NE.0 ) THEN 00502 CALL XERBLA( 'DLATME', -INFO ) 00503 RETURN 00504 END IF 00505 * 00506 * Initialize random number generator 00507 * 00508 DO 30 I = 1, 4 00509 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 ) 00510 30 CONTINUE 00511 * 00512 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 00513 $ ISEED( 4 ) = ISEED( 4 ) + 1 00514 * 00515 * 2) Set up diagonal of A 00516 * 00517 * Compute D according to COND and MODE 00518 * 00519 CALL DLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO ) 00520 IF( IINFO.NE.0 ) THEN 00521 INFO = 1 00522 RETURN 00523 END IF 00524 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN 00525 * 00526 * Scale by DMAX 00527 * 00528 TEMP = ABS( D( 1 ) ) 00529 DO 40 I = 2, N 00530 TEMP = MAX( TEMP, ABS( D( I ) ) ) 00531 40 CONTINUE 00532 * 00533 IF( TEMP.GT.ZERO ) THEN 00534 ALPHA = DMAX / TEMP 00535 ELSE IF( DMAX.NE.ZERO ) THEN 00536 INFO = 2 00537 RETURN 00538 ELSE 00539 ALPHA = ZERO 00540 END IF 00541 * 00542 CALL DSCAL( N, ALPHA, D, 1 ) 00543 * 00544 END IF 00545 * 00546 CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 00547 CALL DCOPY( N, D, 1, A, LDA+1 ) 00548 * 00549 * Set up complex conjugate pairs 00550 * 00551 IF( MODE.EQ.0 ) THEN 00552 IF( USEEI ) THEN 00553 DO 50 J = 2, N 00554 IF( LSAME( EI( J ), 'I' ) ) THEN 00555 A( J-1, J ) = A( J, J ) 00556 A( J, J-1 ) = -A( J, J ) 00557 A( J, J ) = A( J-1, J-1 ) 00558 END IF 00559 50 CONTINUE 00560 END IF 00561 * 00562 ELSE IF( ABS( MODE ).EQ.5 ) THEN 00563 * 00564 DO 60 J = 2, N, 2 00565 IF( DLARAN( ISEED ).GT.HALF ) THEN 00566 A( J-1, J ) = A( J, J ) 00567 A( J, J-1 ) = -A( J, J ) 00568 A( J, J ) = A( J-1, J-1 ) 00569 END IF 00570 60 CONTINUE 00571 END IF 00572 * 00573 * 3) If UPPER='T', set upper triangle of A to random numbers. 00574 * (but don't modify the corners of 2x2 blocks.) 00575 * 00576 IF( IUPPER.NE.0 ) THEN 00577 DO 70 JC = 2, N 00578 IF( A( JC-1, JC ).NE.ZERO ) THEN 00579 JR = JC - 2 00580 ELSE 00581 JR = JC - 1 00582 END IF 00583 CALL DLARNV( IDIST, ISEED, JR, A( 1, JC ) ) 00584 70 CONTINUE 00585 END IF 00586 * 00587 * 4) If SIM='T', apply similarity transformation. 00588 * 00589 * -1 00590 * Transform is X A X , where X = U S V, thus 00591 * 00592 * it is U S V A V' (1/S) U' 00593 * 00594 IF( ISIM.NE.0 ) THEN 00595 * 00596 * Compute S (singular values of the eigenvector matrix) 00597 * according to CONDS and MODES 00598 * 00599 CALL DLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO ) 00600 IF( IINFO.NE.0 ) THEN 00601 INFO = 3 00602 RETURN 00603 END IF 00604 * 00605 * Multiply by V and V' 00606 * 00607 CALL DLARGE( N, A, LDA, ISEED, WORK, IINFO ) 00608 IF( IINFO.NE.0 ) THEN 00609 INFO = 4 00610 RETURN 00611 END IF 00612 * 00613 * Multiply by S and (1/S) 00614 * 00615 DO 80 J = 1, N 00616 CALL DSCAL( N, DS( J ), A( J, 1 ), LDA ) 00617 IF( DS( J ).NE.ZERO ) THEN 00618 CALL DSCAL( N, ONE / DS( J ), A( 1, J ), 1 ) 00619 ELSE 00620 INFO = 5 00621 RETURN 00622 END IF 00623 80 CONTINUE 00624 * 00625 * Multiply by U and U' 00626 * 00627 CALL DLARGE( N, A, LDA, ISEED, WORK, IINFO ) 00628 IF( IINFO.NE.0 ) THEN 00629 INFO = 4 00630 RETURN 00631 END IF 00632 END IF 00633 * 00634 * 5) Reduce the bandwidth. 00635 * 00636 IF( KL.LT.N-1 ) THEN 00637 * 00638 * Reduce bandwidth -- kill column 00639 * 00640 DO 90 JCR = KL + 1, N - 1 00641 IC = JCR - KL 00642 IROWS = N + 1 - JCR 00643 ICOLS = N + KL - JCR 00644 * 00645 CALL DCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 ) 00646 XNORMS = WORK( 1 ) 00647 CALL DLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU ) 00648 WORK( 1 ) = ONE 00649 * 00650 CALL DGEMV( 'T', IROWS, ICOLS, ONE, A( JCR, IC+1 ), LDA, 00651 $ WORK, 1, ZERO, WORK( IROWS+1 ), 1 ) 00652 CALL DGER( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1, 00653 $ A( JCR, IC+1 ), LDA ) 00654 * 00655 CALL DGEMV( 'N', N, IROWS, ONE, A( 1, JCR ), LDA, WORK, 1, 00656 $ ZERO, WORK( IROWS+1 ), 1 ) 00657 CALL DGER( N, IROWS, -TAU, WORK( IROWS+1 ), 1, WORK, 1, 00658 $ A( 1, JCR ), LDA ) 00659 * 00660 A( JCR, IC ) = XNORMS 00661 CALL DLASET( 'Full', IROWS-1, 1, ZERO, ZERO, A( JCR+1, IC ), 00662 $ LDA ) 00663 90 CONTINUE 00664 ELSE IF( KU.LT.N-1 ) THEN 00665 * 00666 * Reduce upper bandwidth -- kill a row at a time. 00667 * 00668 DO 100 JCR = KU + 1, N - 1 00669 IR = JCR - KU 00670 IROWS = N + KU - JCR 00671 ICOLS = N + 1 - JCR 00672 * 00673 CALL DCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 ) 00674 XNORMS = WORK( 1 ) 00675 CALL DLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU ) 00676 WORK( 1 ) = ONE 00677 * 00678 CALL DGEMV( 'N', IROWS, ICOLS, ONE, A( IR+1, JCR ), LDA, 00679 $ WORK, 1, ZERO, WORK( ICOLS+1 ), 1 ) 00680 CALL DGER( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1, 00681 $ A( IR+1, JCR ), LDA ) 00682 * 00683 CALL DGEMV( 'C', ICOLS, N, ONE, A( JCR, 1 ), LDA, WORK, 1, 00684 $ ZERO, WORK( ICOLS+1 ), 1 ) 00685 CALL DGER( ICOLS, N, -TAU, WORK, 1, WORK( ICOLS+1 ), 1, 00686 $ A( JCR, 1 ), LDA ) 00687 * 00688 A( IR, JCR ) = XNORMS 00689 CALL DLASET( 'Full', 1, ICOLS-1, ZERO, ZERO, A( IR, JCR+1 ), 00690 $ LDA ) 00691 100 CONTINUE 00692 END IF 00693 * 00694 * Scale the matrix to have norm ANORM 00695 * 00696 IF( ANORM.GE.ZERO ) THEN 00697 TEMP = DLANGE( 'M', N, N, A, LDA, TEMPA ) 00698 IF( TEMP.GT.ZERO ) THEN 00699 ALPHA = ANORM / TEMP 00700 DO 110 J = 1, N 00701 CALL DSCAL( N, ALPHA, A( 1, J ), 1 ) 00702 110 CONTINUE 00703 END IF 00704 END IF 00705 * 00706 RETURN 00707 * 00708 * End of DLATME 00709 * 00710 END