1*> \brief \b SLATM4 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE SLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, 12* TRIANG, IDIST, ISEED, A, LDA ) 13* 14* .. Scalar Arguments .. 15* INTEGER IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2 16* REAL AMAGN, RCOND, TRIANG 17* .. 18* .. Array Arguments .. 19* INTEGER ISEED( 4 ) 20* REAL A( LDA, * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> SLATM4 generates basic square matrices, which may later be 30*> multiplied by others in order to produce test matrices. It is 31*> intended mainly to be used to test the generalized eigenvalue 32*> routines. 33*> 34*> It first generates the diagonal and (possibly) subdiagonal, 35*> according to the value of ITYPE, NZ1, NZ2, ISIGN, AMAGN, and RCOND. 36*> It then fills in the upper triangle with random numbers, if TRIANG is 37*> non-zero. 38*> \endverbatim 39* 40* Arguments: 41* ========== 42* 43*> \param[in] ITYPE 44*> \verbatim 45*> ITYPE is INTEGER 46*> The "type" of matrix on the diagonal and sub-diagonal. 47*> If ITYPE < 0, then type abs(ITYPE) is generated and then 48*> swapped end for end (A(I,J) := A'(N-J,N-I).) See also 49*> the description of AMAGN and ISIGN. 50*> 51*> Special types: 52*> = 0: the zero matrix. 53*> = 1: the identity. 54*> = 2: a transposed Jordan block. 55*> = 3: If N is odd, then a k+1 x k+1 transposed Jordan block 56*> followed by a k x k identity block, where k=(N-1)/2. 57*> If N is even, then k=(N-2)/2, and a zero diagonal entry 58*> is tacked onto the end. 59*> 60*> Diagonal types. The diagonal consists of NZ1 zeros, then 61*> k=N-NZ1-NZ2 nonzeros. The subdiagonal is zero. ITYPE 62*> specifies the nonzero diagonal entries as follows: 63*> = 4: 1, ..., k 64*> = 5: 1, RCOND, ..., RCOND 65*> = 6: 1, ..., 1, RCOND 66*> = 7: 1, a, a^2, ..., a^(k-1)=RCOND 67*> = 8: 1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND 68*> = 9: random numbers chosen from (RCOND,1) 69*> = 10: random numbers with distribution IDIST (see SLARND.) 70*> \endverbatim 71*> 72*> \param[in] N 73*> \verbatim 74*> N is INTEGER 75*> The order of the matrix. 76*> \endverbatim 77*> 78*> \param[in] NZ1 79*> \verbatim 80*> NZ1 is INTEGER 81*> If abs(ITYPE) > 3, then the first NZ1 diagonal entries will 82*> be zero. 83*> \endverbatim 84*> 85*> \param[in] NZ2 86*> \verbatim 87*> NZ2 is INTEGER 88*> If abs(ITYPE) > 3, then the last NZ2 diagonal entries will 89*> be zero. 90*> \endverbatim 91*> 92*> \param[in] ISIGN 93*> \verbatim 94*> ISIGN is INTEGER 95*> = 0: The sign of the diagonal and subdiagonal entries will 96*> be left unchanged. 97*> = 1: The diagonal and subdiagonal entries will have their 98*> sign changed at random. 99*> = 2: If ITYPE is 2 or 3, then the same as ISIGN=1. 100*> Otherwise, with probability 0.5, odd-even pairs of 101*> diagonal entries A(2*j-1,2*j-1), A(2*j,2*j) will be 102*> converted to a 2x2 block by pre- and post-multiplying 103*> by distinct random orthogonal rotations. The remaining 104*> diagonal entries will have their sign changed at random. 105*> \endverbatim 106*> 107*> \param[in] AMAGN 108*> \verbatim 109*> AMAGN is REAL 110*> The diagonal and subdiagonal entries will be multiplied by 111*> AMAGN. 112*> \endverbatim 113*> 114*> \param[in] RCOND 115*> \verbatim 116*> RCOND is REAL 117*> If abs(ITYPE) > 4, then the smallest diagonal entry will be 118*> entry will be RCOND. RCOND must be between 0 and 1. 119*> \endverbatim 120*> 121*> \param[in] TRIANG 122*> \verbatim 123*> TRIANG is REAL 124*> The entries above the diagonal will be random numbers with 125*> magnitude bounded by TRIANG (i.e., random numbers multiplied 126*> by TRIANG.) 127*> \endverbatim 128*> 129*> \param[in] IDIST 130*> \verbatim 131*> IDIST is INTEGER 132*> Specifies the type of distribution to be used to generate a 133*> random matrix. 134*> = 1: UNIFORM( 0, 1 ) 135*> = 2: UNIFORM( -1, 1 ) 136*> = 3: NORMAL ( 0, 1 ) 137*> \endverbatim 138*> 139*> \param[in,out] ISEED 140*> \verbatim 141*> ISEED is INTEGER array, dimension (4) 142*> On entry ISEED specifies the seed of the random number 143*> generator. The values of ISEED are changed on exit, and can 144*> be used in the next call to SLATM4 to continue the same 145*> random number sequence. 146*> Note: ISEED(4) should be odd, for the random number generator 147*> used at present. 148*> \endverbatim 149*> 150*> \param[out] A 151*> \verbatim 152*> A is REAL array, dimension (LDA, N) 153*> Array to be computed. 154*> \endverbatim 155*> 156*> \param[in] LDA 157*> \verbatim 158*> LDA is INTEGER 159*> Leading dimension of A. Must be at least 1 and at least N. 160*> \endverbatim 161* 162* Authors: 163* ======== 164* 165*> \author Univ. of Tennessee 166*> \author Univ. of California Berkeley 167*> \author Univ. of Colorado Denver 168*> \author NAG Ltd. 169* 170*> \date December 2016 171* 172*> \ingroup single_eig 173* 174* ===================================================================== 175 SUBROUTINE SLATM4( ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, 176 $ TRIANG, IDIST, ISEED, A, LDA ) 177* 178* -- LAPACK test routine (version 3.7.0) -- 179* -- LAPACK is a software package provided by Univ. of Tennessee, -- 180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 181* December 2016 182* 183* .. Scalar Arguments .. 184 INTEGER IDIST, ISIGN, ITYPE, LDA, N, NZ1, NZ2 185 REAL AMAGN, RCOND, TRIANG 186* .. 187* .. Array Arguments .. 188 INTEGER ISEED( 4 ) 189 REAL A( LDA, * ) 190* .. 191* 192* ===================================================================== 193* 194* .. Parameters .. 195 REAL ZERO, ONE, TWO 196 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 ) 197 REAL HALF 198 PARAMETER ( HALF = ONE / TWO ) 199* .. 200* .. Local Scalars .. 201 INTEGER I, IOFF, ISDB, ISDE, JC, JD, JR, K, KBEG, KEND, 202 $ KLEN 203 REAL ALPHA, CL, CR, SAFMIN, SL, SR, SV1, SV2, TEMP 204* .. 205* .. External Functions .. 206 REAL SLAMCH, SLARAN, SLARND 207 EXTERNAL SLAMCH, SLARAN, SLARND 208* .. 209* .. External Subroutines .. 210 EXTERNAL SLASET 211* .. 212* .. Intrinsic Functions .. 213 INTRINSIC ABS, EXP, LOG, MAX, MIN, MOD, REAL, SQRT 214* .. 215* .. Executable Statements .. 216* 217 IF( N.LE.0 ) 218 $ RETURN 219 CALL SLASET( 'Full', N, N, ZERO, ZERO, A, LDA ) 220* 221* Insure a correct ISEED 222* 223 IF( MOD( ISEED( 4 ), 2 ).NE.1 ) 224 $ ISEED( 4 ) = ISEED( 4 ) + 1 225* 226* Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2, 227* and RCOND 228* 229 IF( ITYPE.NE.0 ) THEN 230 IF( ABS( ITYPE ).GE.4 ) THEN 231 KBEG = MAX( 1, MIN( N, NZ1+1 ) ) 232 KEND = MAX( KBEG, MIN( N, N-NZ2 ) ) 233 KLEN = KEND + 1 - KBEG 234 ELSE 235 KBEG = 1 236 KEND = N 237 KLEN = N 238 END IF 239 ISDB = 1 240 ISDE = 0 241 GO TO ( 10, 30, 50, 80, 100, 120, 140, 160, 242 $ 180, 200 )ABS( ITYPE ) 243* 244* abs(ITYPE) = 1: Identity 245* 246 10 CONTINUE 247 DO 20 JD = 1, N 248 A( JD, JD ) = ONE 249 20 CONTINUE 250 GO TO 220 251* 252* abs(ITYPE) = 2: Transposed Jordan block 253* 254 30 CONTINUE 255 DO 40 JD = 1, N - 1 256 A( JD+1, JD ) = ONE 257 40 CONTINUE 258 ISDB = 1 259 ISDE = N - 1 260 GO TO 220 261* 262* abs(ITYPE) = 3: Transposed Jordan block, followed by the 263* identity. 264* 265 50 CONTINUE 266 K = ( N-1 ) / 2 267 DO 60 JD = 1, K 268 A( JD+1, JD ) = ONE 269 60 CONTINUE 270 ISDB = 1 271 ISDE = K 272 DO 70 JD = K + 2, 2*K + 1 273 A( JD, JD ) = ONE 274 70 CONTINUE 275 GO TO 220 276* 277* abs(ITYPE) = 4: 1,...,k 278* 279 80 CONTINUE 280 DO 90 JD = KBEG, KEND 281 A( JD, JD ) = REAL( JD-NZ1 ) 282 90 CONTINUE 283 GO TO 220 284* 285* abs(ITYPE) = 5: One large D value: 286* 287 100 CONTINUE 288 DO 110 JD = KBEG + 1, KEND 289 A( JD, JD ) = RCOND 290 110 CONTINUE 291 A( KBEG, KBEG ) = ONE 292 GO TO 220 293* 294* abs(ITYPE) = 6: One small D value: 295* 296 120 CONTINUE 297 DO 130 JD = KBEG, KEND - 1 298 A( JD, JD ) = ONE 299 130 CONTINUE 300 A( KEND, KEND ) = RCOND 301 GO TO 220 302* 303* abs(ITYPE) = 7: Exponentially distributed D values: 304* 305 140 CONTINUE 306 A( KBEG, KBEG ) = ONE 307 IF( KLEN.GT.1 ) THEN 308 ALPHA = RCOND**( ONE / REAL( KLEN-1 ) ) 309 DO 150 I = 2, KLEN 310 A( NZ1+I, NZ1+I ) = ALPHA**REAL( I-1 ) 311 150 CONTINUE 312 END IF 313 GO TO 220 314* 315* abs(ITYPE) = 8: Arithmetically distributed D values: 316* 317 160 CONTINUE 318 A( KBEG, KBEG ) = ONE 319 IF( KLEN.GT.1 ) THEN 320 ALPHA = ( ONE-RCOND ) / REAL( KLEN-1 ) 321 DO 170 I = 2, KLEN 322 A( NZ1+I, NZ1+I ) = REAL( KLEN-I )*ALPHA + RCOND 323 170 CONTINUE 324 END IF 325 GO TO 220 326* 327* abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1): 328* 329 180 CONTINUE 330 ALPHA = LOG( RCOND ) 331 DO 190 JD = KBEG, KEND 332 A( JD, JD ) = EXP( ALPHA*SLARAN( ISEED ) ) 333 190 CONTINUE 334 GO TO 220 335* 336* abs(ITYPE) = 10: Randomly distributed D values from DIST 337* 338 200 CONTINUE 339 DO 210 JD = KBEG, KEND 340 A( JD, JD ) = SLARND( IDIST, ISEED ) 341 210 CONTINUE 342* 343 220 CONTINUE 344* 345* Scale by AMAGN 346* 347 DO 230 JD = KBEG, KEND 348 A( JD, JD ) = AMAGN*REAL( A( JD, JD ) ) 349 230 CONTINUE 350 DO 240 JD = ISDB, ISDE 351 A( JD+1, JD ) = AMAGN*REAL( A( JD+1, JD ) ) 352 240 CONTINUE 353* 354* If ISIGN = 1 or 2, assign random signs to diagonal and 355* subdiagonal 356* 357 IF( ISIGN.GT.0 ) THEN 358 DO 250 JD = KBEG, KEND 359 IF( REAL( A( JD, JD ) ).NE.ZERO ) THEN 360 IF( SLARAN( ISEED ).GT.HALF ) 361 $ A( JD, JD ) = -A( JD, JD ) 362 END IF 363 250 CONTINUE 364 DO 260 JD = ISDB, ISDE 365 IF( REAL( A( JD+1, JD ) ).NE.ZERO ) THEN 366 IF( SLARAN( ISEED ).GT.HALF ) 367 $ A( JD+1, JD ) = -A( JD+1, JD ) 368 END IF 369 260 CONTINUE 370 END IF 371* 372* Reverse if ITYPE < 0 373* 374 IF( ITYPE.LT.0 ) THEN 375 DO 270 JD = KBEG, ( KBEG+KEND-1 ) / 2 376 TEMP = A( JD, JD ) 377 A( JD, JD ) = A( KBEG+KEND-JD, KBEG+KEND-JD ) 378 A( KBEG+KEND-JD, KBEG+KEND-JD ) = TEMP 379 270 CONTINUE 380 DO 280 JD = 1, ( N-1 ) / 2 381 TEMP = A( JD+1, JD ) 382 A( JD+1, JD ) = A( N+1-JD, N-JD ) 383 A( N+1-JD, N-JD ) = TEMP 384 280 CONTINUE 385 END IF 386* 387* If ISIGN = 2, and no subdiagonals already, then apply 388* random rotations to make 2x2 blocks. 389* 390 IF( ISIGN.EQ.2 .AND. ITYPE.NE.2 .AND. ITYPE.NE.3 ) THEN 391 SAFMIN = SLAMCH( 'S' ) 392 DO 290 JD = KBEG, KEND - 1, 2 393 IF( SLARAN( ISEED ).GT.HALF ) THEN 394* 395* Rotation on left. 396* 397 CL = TWO*SLARAN( ISEED ) - ONE 398 SL = TWO*SLARAN( ISEED ) - ONE 399 TEMP = ONE / MAX( SAFMIN, SQRT( CL**2+SL**2 ) ) 400 CL = CL*TEMP 401 SL = SL*TEMP 402* 403* Rotation on right. 404* 405 CR = TWO*SLARAN( ISEED ) - ONE 406 SR = TWO*SLARAN( ISEED ) - ONE 407 TEMP = ONE / MAX( SAFMIN, SQRT( CR**2+SR**2 ) ) 408 CR = CR*TEMP 409 SR = SR*TEMP 410* 411* Apply 412* 413 SV1 = A( JD, JD ) 414 SV2 = A( JD+1, JD+1 ) 415 A( JD, JD ) = CL*CR*SV1 + SL*SR*SV2 416 A( JD+1, JD ) = -SL*CR*SV1 + CL*SR*SV2 417 A( JD, JD+1 ) = -CL*SR*SV1 + SL*CR*SV2 418 A( JD+1, JD+1 ) = SL*SR*SV1 + CL*CR*SV2 419 END IF 420 290 CONTINUE 421 END IF 422* 423 END IF 424* 425* Fill in upper triangle (except for 2x2 blocks) 426* 427 IF( TRIANG.NE.ZERO ) THEN 428 IF( ISIGN.NE.2 .OR. ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN 429 IOFF = 1 430 ELSE 431 IOFF = 2 432 DO 300 JR = 1, N - 1 433 IF( A( JR+1, JR ).EQ.ZERO ) 434 $ A( JR, JR+1 ) = TRIANG*SLARND( IDIST, ISEED ) 435 300 CONTINUE 436 END IF 437* 438 DO 320 JC = 2, N 439 DO 310 JR = 1, JC - IOFF 440 A( JR, JC ) = TRIANG*SLARND( IDIST, ISEED ) 441 310 CONTINUE 442 320 CONTINUE 443 END IF 444* 445 RETURN 446* 447* End of SLATM4 448* 449 END 450