![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLATM4 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 DLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, 00012 * TRIANG, IDIST, ISEED, A, LDA ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2 00016 * DOUBLE PRECISION AMAGN, RCOND, TRIANG 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER ISEED( 4 ) 00020 * DOUBLE PRECISION A( LDA, * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> DLATM4 generates basic square matrices, which may later be 00030 *> multiplied by others in order to produce test matrices. It is 00031 *> intended mainly to be used to test the generalized eigenvalue 00032 *> routines. 00033 *> 00034 *> It first generates the diagonal and (possibly) subdiagonal, 00035 *> according to the value of ITYPE, NZ1, NZ2, ISIGN, AMAGN, and RCOND. 00036 *> It then fills in the upper triangle with random numbers, if TRIANG is 00037 *> non-zero. 00038 *> \endverbatim 00039 * 00040 * Arguments: 00041 * ========== 00042 * 00043 *> \param[in] ITYPE 00044 *> \verbatim 00045 *> ITYPE is INTEGER 00046 *> The "type" of matrix on the diagonal and sub-diagonal. 00047 *> If ITYPE < 0, then type abs(ITYPE) is generated and then 00048 *> swapped end for end (A(I,J) := A'(N-J,N-I).) See also 00049 *> the description of AMAGN and ISIGN. 00050 *> 00051 *> Special types: 00052 *> = 0: the zero matrix. 00053 *> = 1: the identity. 00054 *> = 2: a transposed Jordan block. 00055 *> = 3: If N is odd, then a k+1 x k+1 transposed Jordan block 00056 *> followed by a k x k identity block, where k=(N-1)/2. 00057 *> If N is even, then k=(N-2)/2, and a zero diagonal entry 00058 *> is tacked onto the end. 00059 *> 00060 *> Diagonal types. The diagonal consists of NZ1 zeros, then 00061 *> k=N-NZ1-NZ2 nonzeros. The subdiagonal is zero. ITYPE 00062 *> specifies the nonzero diagonal entries as follows: 00063 *> = 4: 1, ..., k 00064 *> = 5: 1, RCOND, ..., RCOND 00065 *> = 6: 1, ..., 1, RCOND 00066 *> = 7: 1, a, a^2, ..., a^(k-1)=RCOND 00067 *> = 8: 1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND 00068 *> = 9: random numbers chosen from (RCOND,1) 00069 *> = 10: random numbers with distribution IDIST (see DLARND.) 00070 *> \endverbatim 00071 *> 00072 *> \param[in] N 00073 *> \verbatim 00074 *> N is INTEGER 00075 *> The order of the matrix. 00076 *> \endverbatim 00077 *> 00078 *> \param[in] NZ1 00079 *> \verbatim 00080 *> NZ1 is INTEGER 00081 *> If abs(ITYPE) > 3, then the first NZ1 diagonal entries will 00082 *> be zero. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] NZ2 00086 *> \verbatim 00087 *> NZ2 is INTEGER 00088 *> If abs(ITYPE) > 3, then the last NZ2 diagonal entries will 00089 *> be zero. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] ISIGN 00093 *> \verbatim 00094 *> ISIGN is INTEGER 00095 *> = 0: The sign of the diagonal and subdiagonal entries will 00096 *> be left unchanged. 00097 *> = 1: The diagonal and subdiagonal entries will have their 00098 *> sign changed at random. 00099 *> = 2: If ITYPE is 2 or 3, then the same as ISIGN=1. 00100 *> Otherwise, with probability 0.5, odd-even pairs of 00101 *> diagonal entries A(2*j-1,2*j-1), A(2*j,2*j) will be 00102 *> converted to a 2x2 block by pre- and post-multiplying 00103 *> by distinct random orthogonal rotations. The remaining 00104 *> diagonal entries will have their sign changed at random. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] AMAGN 00108 *> \verbatim 00109 *> AMAGN is DOUBLE PRECISION 00110 *> The diagonal and subdiagonal entries will be multiplied by 00111 *> AMAGN. 00112 *> \endverbatim 00113 *> 00114 *> \param[in] RCOND 00115 *> \verbatim 00116 *> RCOND is DOUBLE PRECISION 00117 *> If abs(ITYPE) > 4, then the smallest diagonal entry will be 00118 *> entry will be RCOND. RCOND must be between 0 and 1. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] TRIANG 00122 *> \verbatim 00123 *> TRIANG is DOUBLE PRECISION 00124 *> The entries above the diagonal will be random numbers with 00125 *> magnitude bounded by TRIANG (i.e., random numbers multiplied 00126 *> by TRIANG.) 00127 *> \endverbatim 00128 *> 00129 *> \param[in] IDIST 00130 *> \verbatim 00131 *> IDIST is INTEGER 00132 *> Specifies the type of distribution to be used to generate a 00133 *> random matrix. 00134 *> = 1: UNIFORM( 0, 1 ) 00135 *> = 2: UNIFORM( -1, 1 ) 00136 *> = 3: NORMAL ( 0, 1 ) 00137 *> \endverbatim 00138 *> 00139 *> \param[in,out] ISEED 00140 *> \verbatim 00141 *> ISEED is INTEGER array, dimension (4) 00142 *> On entry ISEED specifies the seed of the random number 00143 *> generator. The values of ISEED are changed on exit, and can 00144 *> be used in the next call to DLATM4 to continue the same 00145 *> random number sequence. 00146 *> Note: ISEED(4) should be odd, for the random number generator 00147 *> used at present. 00148 *> \endverbatim 00149 *> 00150 *> \param[out] A 00151 *> \verbatim 00152 *> A is DOUBLE PRECISION array, dimension (LDA, N) 00153 *> Array to be computed. 00154 *> \endverbatim 00155 *> 00156 *> \param[in] LDA 00157 *> \verbatim 00158 *> LDA is INTEGER 00159 *> Leading dimension of A. Must be at least 1 and at least N. 00160 *> \endverbatim 00161 * 00162 * Authors: 00163 * ======== 00164 * 00165 *> \author Univ. of Tennessee 00166 *> \author Univ. of California Berkeley 00167 *> \author Univ. of Colorado Denver 00168 *> \author NAG Ltd. 00169 * 00170 *> \date November 2011 00171 * 00172 *> \ingroup double_eig 00173 * 00174 * ===================================================================== 00175 SUBROUTINE DLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, 00176 $ TRIANG, IDIST, ISEED, A, LDA ) 00177 * 00178 * -- LAPACK test routine (version 3.4.0) -- 00179 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00180 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00181 * November 2011 00182 * 00183 * .. Scalar Arguments .. 00184 INTEGER IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2 00185 DOUBLE PRECISION AMAGN, RCOND, TRIANG 00186 * .. 00187 * .. Array Arguments .. 00188 INTEGER ISEED( 4 ) 00189 DOUBLE PRECISION A( LDA, * ) 00190 * .. 00191 * 00192 * ===================================================================== 00193 * 00194 * .. Parameters .. 00195 DOUBLE PRECISION ZERO, ONE, TWO 00196 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00197 DOUBLE PRECISION HALF 00198 PARAMETER ( HALF = ONE / TWO ) 00199 * .. 00200 * .. Local Scalars .. 00201 INTEGER I, IOFF, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND, 00202 $ KLEN 00203 DOUBLE PRECISION ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP 00204 * .. 00205 * .. External Functions .. 00206 DOUBLE PRECISION DLAMCH, DLARAN, DLARND 00207 EXTERNAL DLAMCH, DLARAN, DLARND 00208 * .. 00209 * .. External Subroutines .. 00210 EXTERNAL DLASET 00211 * .. 00212 * .. Intrinsic Functions .. 00213 INTRINSIC ABS, DBLE, EXP, LOG, MAX, MIN, MOD, SQRT 00214 * .. 00215 * .. Executable Statements .. 00216 * 00217 IF( N.LE.0 ) 00218 $ RETURN 00219 CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 00220 * 00221 * Insure a correct ISEED 00222 * 00223 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 00224 $ ISEED( 4 ) = ISEED( 4 ) + 1 00225 * 00226 * Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2, 00227 * and RCOND 00228 * 00229 IF( ITYPE.NE.0 ) THEN 00230 IF( ABS( ITYPE ).GE.4 ) THEN 00231 KBEG = MAX( 1, MIN( N, NZ1+1 ) ) 00232 KEND = MAX( KBEG, MIN( N, N-NZ2 ) ) 00233 KLEN = KEND + 1 - KBEG 00234 ELSE 00235 KBEG = 1 00236 KEND = N 00237 KLEN = N 00238 END IF 00239 ISDB = 1 00240 ISDE = 0 00241 GO TO ( 10, 30, 50, 80, 100, 120, 140, 160, 00242 $ 180, 200 )ABS( ITYPE ) 00243 * 00244 * abs(ITYPE) = 1: Identity 00245 * 00246 10 CONTINUE 00247 DO 20 JD = 1, N 00248 A( JD, JD ) = ONE 00249 20 CONTINUE 00250 GO TO 220 00251 * 00252 * abs(ITYPE) = 2: Transposed Jordan block 00253 * 00254 30 CONTINUE 00255 DO 40 JD = 1, N - 1 00256 A( JD+1, JD ) = ONE 00257 40 CONTINUE 00258 ISDB = 1 00259 ISDE = N - 1 00260 GO TO 220 00261 * 00262 * abs(ITYPE) = 3: Transposed Jordan block, followed by the 00263 * identity. 00264 * 00265 50 CONTINUE 00266 K = ( N-1 ) / 2 00267 DO 60 JD = 1, K 00268 A( JD+1, JD ) = ONE 00269 60 CONTINUE 00270 ISDB = 1 00271 ISDE = K 00272 DO 70 JD = K + 2, 2*K + 1 00273 A( JD, JD ) = ONE 00274 70 CONTINUE 00275 GO TO 220 00276 * 00277 * abs(ITYPE) = 4: 1,...,k 00278 * 00279 80 CONTINUE 00280 DO 90 JD = KBEG, KEND 00281 A( JD, JD ) = DBLE( JD-NZ1 ) 00282 90 CONTINUE 00283 GO TO 220 00284 * 00285 * abs(ITYPE) = 5: One large D value: 00286 * 00287 100 CONTINUE 00288 DO 110 JD = KBEG + 1, KEND 00289 A( JD, JD ) = RCOND 00290 110 CONTINUE 00291 A( KBEG, KBEG ) = ONE 00292 GO TO 220 00293 * 00294 * abs(ITYPE) = 6: One small D value: 00295 * 00296 120 CONTINUE 00297 DO 130 JD = KBEG, KEND - 1 00298 A( JD, JD ) = ONE 00299 130 CONTINUE 00300 A( KEND, KEND ) = RCOND 00301 GO TO 220 00302 * 00303 * abs(ITYPE) = 7: Exponentially distributed D values: 00304 * 00305 140 CONTINUE 00306 A( KBEG, KBEG ) = ONE 00307 IF( KLEN.GT.1 ) THEN 00308 ALPHA = RCOND**( ONE / DBLE( KLEN-1 ) ) 00309 DO 150 I = 2, KLEN 00310 A( NZ1+I, NZ1+I ) = ALPHA**DBLE( I-1 ) 00311 150 CONTINUE 00312 END IF 00313 GO TO 220 00314 * 00315 * abs(ITYPE) = 8: Arithmetically distributed D values: 00316 * 00317 160 CONTINUE 00318 A( KBEG, KBEG ) = ONE 00319 IF( KLEN.GT.1 ) THEN 00320 ALPHA = ( ONE-RCOND ) / DBLE( KLEN-1 ) 00321 DO 170 I = 2, KLEN 00322 A( NZ1+I, NZ1+I ) = DBLE( KLEN-I )*ALPHA + RCOND 00323 170 CONTINUE 00324 END IF 00325 GO TO 220 00326 * 00327 * abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1): 00328 * 00329 180 CONTINUE 00330 ALPHA = LOG( RCOND ) 00331 DO 190 JD = KBEG, KEND 00332 A( JD, JD ) = EXP( ALPHA*DLARAN( ISEED ) ) 00333 190 CONTINUE 00334 GO TO 220 00335 * 00336 * abs(ITYPE) = 10: Randomly distributed D values from DIST 00337 * 00338 200 CONTINUE 00339 DO 210 JD = KBEG, KEND 00340 A( JD, JD ) = DLARND( IDIST, ISEED ) 00341 210 CONTINUE 00342 * 00343 220 CONTINUE 00344 * 00345 * Scale by AMAGN 00346 * 00347 DO 230 JD = KBEG, KEND 00348 A( JD, JD ) = AMAGN*DBLE( A( JD, JD ) ) 00349 230 CONTINUE 00350 DO 240 JD = ISDB, ISDE 00351 A( JD+1, JD ) = AMAGN*DBLE( A( JD+1, JD ) ) 00352 240 CONTINUE 00353 * 00354 * If ISIGN = 1 or 2, assign random signs to diagonal and 00355 * subdiagonal 00356 * 00357 IF( ISIGN.GT.0 ) THEN 00358 DO 250 JD = KBEG, KEND 00359 IF( DBLE( A( JD, JD ) ).NE.ZERO ) THEN 00360 IF( DLARAN( ISEED ).GT.HALF ) 00361 $ A( JD, JD ) = -A( JD, JD ) 00362 END IF 00363 250 CONTINUE 00364 DO 260 JD = ISDB, ISDE 00365 IF( DBLE( A( JD+1, JD ) ).NE.ZERO ) THEN 00366 IF( DLARAN( ISEED ).GT.HALF ) 00367 $ A( JD+1, JD ) = -A( JD+1, JD ) 00368 END IF 00369 260 CONTINUE 00370 END IF 00371 * 00372 * Reverse if ITYPE < 0 00373 * 00374 IF( ITYPE.LT.0 ) THEN 00375 DO 270 JD = KBEG, ( KBEG+KEND-1 ) / 2 00376 TEMP = A( JD, JD ) 00377 A( JD, JD ) = A( KBEG+KEND-JD, KBEG+KEND-JD ) 00378 A( KBEG+KEND-JD, KBEG+KEND-JD ) = TEMP 00379 270 CONTINUE 00380 DO 280 JD = 1, ( N-1 ) / 2 00381 TEMP = A( JD+1, JD ) 00382 A( JD+1, JD ) = A( N+1-JD, N-JD ) 00383 A( N+1-JD, N-JD ) = TEMP 00384 280 CONTINUE 00385 END IF 00386 * 00387 * If ISIGN = 2, and no subdiagonals already, then apply 00388 * random rotations to make 2x2 blocks. 00389 * 00390 IF( ISIGN.EQ.2 .AND. ITYPE.NE.2 .AND. ITYPE.NE.3 ) THEN 00391 SAFMIN = DLAMCH( 'S' ) 00392 DO 290 JD = KBEG, KEND - 1, 2 00393 IF( DLARAN( ISEED ).GT.HALF ) THEN 00394 * 00395 * Rotation on left. 00396 * 00397 CL = TWO*DLARAN( ISEED ) - ONE 00398 SL = TWO*DLARAN( ISEED ) - ONE 00399 TEMP = ONE / MAX( SAFMIN, SQRT( CL**2+SL**2 ) ) 00400 CL = CL*TEMP 00401 SL = SL*TEMP 00402 * 00403 * Rotation on right. 00404 * 00405 CR = TWO*DLARAN( ISEED ) - ONE 00406 SR = TWO*DLARAN( ISEED ) - ONE 00407 TEMP = ONE / MAX( SAFMIN, SQRT( CR**2+SR**2 ) ) 00408 CR = CR*TEMP 00409 SR = SR*TEMP 00410 * 00411 * Apply 00412 * 00413 SV1 = A( JD, JD ) 00414 SV2 = A( JD+1, JD+1 ) 00415 A( JD, JD ) = CL*CR*SV1 + SL*SR*SV2 00416 A( JD+1, JD ) = -SL*CR*SV1 + CL*SR*SV2 00417 A( JD, JD+1 ) = -CL*SR*SV1 + SL*CR*SV2 00418 A( JD+1, JD+1 ) = SL*SR*SV1 + CL*CR*SV2 00419 END IF 00420 290 CONTINUE 00421 END IF 00422 * 00423 END IF 00424 * 00425 * Fill in upper triangle (except for 2x2 blocks) 00426 * 00427 IF( TRIANG.NE.ZERO ) THEN 00428 IF( ISIGN.NE.2 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN 00429 IOFF = 1 00430 ELSE 00431 IOFF = 2 00432 DO 300 JR = 1, N - 1 00433 IF( A( JR+1, JR ).EQ.ZERO ) 00434 $ A( JR, JR+1 ) = TRIANG*DLARND( IDIST, ISEED ) 00435 300 CONTINUE 00436 END IF 00437 * 00438 DO 320 JC = 2, N 00439 DO 310 JR = 1, JC - IOFF 00440 A( JR, JC ) = TRIANG*DLARND( IDIST, ISEED ) 00441 310 CONTINUE 00442 320 CONTINUE 00443 END IF 00444 * 00445 RETURN 00446 * 00447 * End of DLATM4 00448 * 00449 END