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