|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLATTP( IMAT, UPLO, TRANS, DIAG, ISEED, N, A, B, WORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER DIAG, TRANS, UPLO 00010 INTEGER IMAT, INFO, N 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER ISEED( 4 ) 00014 DOUBLE PRECISION A( * ), B( * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLATTP generates a triangular test matrix in packed storage. 00021 * IMAT and UPLO uniquely specify the properties of the test 00022 * matrix, which is returned in the array AP. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * IMAT (input) INTEGER 00028 * An integer key describing which matrix to generate for this 00029 * path. 00030 * 00031 * UPLO (input) CHARACTER*1 00032 * Specifies whether the matrix A will be upper or lower 00033 * triangular. 00034 * = 'U': Upper triangular 00035 * = 'L': Lower triangular 00036 * 00037 * TRANS (input) CHARACTER*1 00038 * Specifies whether the matrix or its transpose will be used. 00039 * = 'N': No transpose 00040 * = 'T': Transpose 00041 * = 'C': Conjugate transpose (= Transpose) 00042 * 00043 * DIAG (output) CHARACTER*1 00044 * Specifies whether or not the matrix A is unit triangular. 00045 * = 'N': Non-unit triangular 00046 * = 'U': Unit triangular 00047 * 00048 * ISEED (input/output) INTEGER array, dimension (4) 00049 * The seed vector for the random number generator (used in 00050 * DLATMS). Modified on exit. 00051 * 00052 * N (input) INTEGER 00053 * The order of the matrix to be generated. 00054 * 00055 * A (output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 00056 * The upper or lower triangular matrix A, packed columnwise in 00057 * a linear array. The j-th column of A is stored in the array 00058 * AP as follows: 00059 * if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j; 00060 * if UPLO = 'L', 00061 * AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n. 00062 * 00063 * B (output) DOUBLE PRECISION array, dimension (N) 00064 * The right hand side vector, if IMAT > 10. 00065 * 00066 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 00067 * 00068 * INFO (output) INTEGER 00069 * = 0: successful exit 00070 * < 0: if INFO = -k, the k-th argument had an illegal value 00071 * 00072 * ===================================================================== 00073 * 00074 * .. Parameters .. 00075 DOUBLE PRECISION ONE, TWO, ZERO 00076 PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0, ZERO = 0.0D+0 ) 00077 * .. 00078 * .. Local Scalars .. 00079 LOGICAL UPPER 00080 CHARACTER DIST, PACKIT, TYPE 00081 CHARACTER*3 PATH 00082 INTEGER I, IY, J, JC, JCNEXT, JCOUNT, JJ, JL, JR, JX, 00083 $ KL, KU, MODE 00084 DOUBLE PRECISION ANORM, BIGNUM, BNORM, BSCAL, C, CNDNUM, PLUS1, 00085 $ PLUS2, RA, RB, REXP, S, SFAC, SMLNUM, STAR1, 00086 $ STEMP, T, TEXP, TLEFT, TSCAL, ULP, UNFL, X, Y, 00087 $ Z 00088 * .. 00089 * .. External Functions .. 00090 LOGICAL LSAME 00091 INTEGER IDAMAX 00092 DOUBLE PRECISION DLAMCH, DLARND 00093 EXTERNAL LSAME, IDAMAX, DLAMCH, DLARND 00094 * .. 00095 * .. External Subroutines .. 00096 EXTERNAL DLABAD, DLARNV, DLATB4, DLATMS, DROT, DROTG, 00097 $ DSCAL 00098 * .. 00099 * .. Intrinsic Functions .. 00100 INTRINSIC ABS, DBLE, MAX, SIGN, SQRT 00101 * .. 00102 * .. Executable Statements .. 00103 * 00104 PATH( 1: 1 ) = 'Double precision' 00105 PATH( 2: 3 ) = 'TP' 00106 UNFL = DLAMCH( 'Safe minimum' ) 00107 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00108 SMLNUM = UNFL 00109 BIGNUM = ( ONE-ULP ) / SMLNUM 00110 CALL DLABAD( SMLNUM, BIGNUM ) 00111 IF( ( IMAT.GE.7 .AND. IMAT.LE.10 ) .OR. IMAT.EQ.18 ) THEN 00112 DIAG = 'U' 00113 ELSE 00114 DIAG = 'N' 00115 END IF 00116 INFO = 0 00117 * 00118 * Quick return if N.LE.0. 00119 * 00120 IF( N.LE.0 ) 00121 $ RETURN 00122 * 00123 * Call DLATB4 to set parameters for SLATMS. 00124 * 00125 UPPER = LSAME( UPLO, 'U' ) 00126 IF( UPPER ) THEN 00127 CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00128 $ CNDNUM, DIST ) 00129 PACKIT = 'C' 00130 ELSE 00131 CALL DLATB4( PATH, -IMAT, N, N, TYPE, KL, KU, ANORM, MODE, 00132 $ CNDNUM, DIST ) 00133 PACKIT = 'R' 00134 END IF 00135 * 00136 * IMAT <= 6: Non-unit triangular matrix 00137 * 00138 IF( IMAT.LE.6 ) THEN 00139 CALL DLATMS( N, N, DIST, ISEED, TYPE, B, MODE, CNDNUM, ANORM, 00140 $ KL, KU, PACKIT, A, N, WORK, INFO ) 00141 * 00142 * IMAT > 6: Unit triangular matrix 00143 * The diagonal is deliberately set to something other than 1. 00144 * 00145 * IMAT = 7: Matrix is the identity 00146 * 00147 ELSE IF( IMAT.EQ.7 ) THEN 00148 IF( UPPER ) THEN 00149 JC = 1 00150 DO 20 J = 1, N 00151 DO 10 I = 1, J - 1 00152 A( JC+I-1 ) = ZERO 00153 10 CONTINUE 00154 A( JC+J-1 ) = J 00155 JC = JC + J 00156 20 CONTINUE 00157 ELSE 00158 JC = 1 00159 DO 40 J = 1, N 00160 A( JC ) = J 00161 DO 30 I = J + 1, N 00162 A( JC+I-J ) = ZERO 00163 30 CONTINUE 00164 JC = JC + N - J + 1 00165 40 CONTINUE 00166 END IF 00167 * 00168 * IMAT > 7: Non-trivial unit triangular matrix 00169 * 00170 * Generate a unit triangular matrix T with condition CNDNUM by 00171 * forming a triangular matrix with known singular values and 00172 * filling in the zero entries with Givens rotations. 00173 * 00174 ELSE IF( IMAT.LE.10 ) THEN 00175 IF( UPPER ) THEN 00176 JC = 0 00177 DO 60 J = 1, N 00178 DO 50 I = 1, J - 1 00179 A( JC+I ) = ZERO 00180 50 CONTINUE 00181 A( JC+J ) = J 00182 JC = JC + J 00183 60 CONTINUE 00184 ELSE 00185 JC = 1 00186 DO 80 J = 1, N 00187 A( JC ) = J 00188 DO 70 I = J + 1, N 00189 A( JC+I-J ) = ZERO 00190 70 CONTINUE 00191 JC = JC + N - J + 1 00192 80 CONTINUE 00193 END IF 00194 * 00195 * Since the trace of a unit triangular matrix is 1, the product 00196 * of its singular values must be 1. Let s = sqrt(CNDNUM), 00197 * x = sqrt(s) - 1/sqrt(s), y = sqrt(2/(n-2))*x, and z = x**2. 00198 * The following triangular matrix has singular values s, 1, 1, 00199 * ..., 1, 1/s: 00200 * 00201 * 1 y y y ... y y z 00202 * 1 0 0 ... 0 0 y 00203 * 1 0 ... 0 0 y 00204 * . ... . . . 00205 * . . . . 00206 * 1 0 y 00207 * 1 y 00208 * 1 00209 * 00210 * To fill in the zeros, we first multiply by a matrix with small 00211 * condition number of the form 00212 * 00213 * 1 0 0 0 0 ... 00214 * 1 + * 0 0 ... 00215 * 1 + 0 0 0 00216 * 1 + * 0 0 00217 * 1 + 0 0 00218 * ... 00219 * 1 + 0 00220 * 1 0 00221 * 1 00222 * 00223 * Each element marked with a '*' is formed by taking the product 00224 * of the adjacent elements marked with '+'. The '*'s can be 00225 * chosen freely, and the '+'s are chosen so that the inverse of 00226 * T will have elements of the same magnitude as T. If the *'s in 00227 * both T and inv(T) have small magnitude, T is well conditioned. 00228 * The two offdiagonals of T are stored in WORK. 00229 * 00230 * The product of these two matrices has the form 00231 * 00232 * 1 y y y y y . y y z 00233 * 1 + * 0 0 . 0 0 y 00234 * 1 + 0 0 . 0 0 y 00235 * 1 + * . . . . 00236 * 1 + . . . . 00237 * . . . . . 00238 * . . . . 00239 * 1 + y 00240 * 1 y 00241 * 1 00242 * 00243 * Now we multiply by Givens rotations, using the fact that 00244 * 00245 * [ c s ] [ 1 w ] [ -c -s ] = [ 1 -w ] 00246 * [ -s c ] [ 0 1 ] [ s -c ] [ 0 1 ] 00247 * and 00248 * [ -c -s ] [ 1 0 ] [ c s ] = [ 1 0 ] 00249 * [ s -c ] [ w 1 ] [ -s c ] [ -w 1 ] 00250 * 00251 * where c = w / sqrt(w**2+4) and s = 2 / sqrt(w**2+4). 00252 * 00253 STAR1 = 0.25D0 00254 SFAC = 0.5D0 00255 PLUS1 = SFAC 00256 DO 90 J = 1, N, 2 00257 PLUS2 = STAR1 / PLUS1 00258 WORK( J ) = PLUS1 00259 WORK( N+J ) = STAR1 00260 IF( J+1.LE.N ) THEN 00261 WORK( J+1 ) = PLUS2 00262 WORK( N+J+1 ) = ZERO 00263 PLUS1 = STAR1 / PLUS2 00264 REXP = DLARND( 2, ISEED ) 00265 STAR1 = STAR1*( SFAC**REXP ) 00266 IF( REXP.LT.ZERO ) THEN 00267 STAR1 = -SFAC**( ONE-REXP ) 00268 ELSE 00269 STAR1 = SFAC**( ONE+REXP ) 00270 END IF 00271 END IF 00272 90 CONTINUE 00273 * 00274 X = SQRT( CNDNUM ) - ONE / SQRT( CNDNUM ) 00275 IF( N.GT.2 ) THEN 00276 Y = SQRT( TWO / DBLE( N-2 ) )*X 00277 ELSE 00278 Y = ZERO 00279 END IF 00280 Z = X*X 00281 * 00282 IF( UPPER ) THEN 00283 * 00284 * Set the upper triangle of A with a unit triangular matrix 00285 * of known condition number. 00286 * 00287 JC = 1 00288 DO 100 J = 2, N 00289 A( JC+1 ) = Y 00290 IF( J.GT.2 ) 00291 $ A( JC+J-1 ) = WORK( J-2 ) 00292 IF( J.GT.3 ) 00293 $ A( JC+J-2 ) = WORK( N+J-3 ) 00294 JC = JC + J 00295 100 CONTINUE 00296 JC = JC - N 00297 A( JC+1 ) = Z 00298 DO 110 J = 2, N - 1 00299 A( JC+J ) = Y 00300 110 CONTINUE 00301 ELSE 00302 * 00303 * Set the lower triangle of A with a unit triangular matrix 00304 * of known condition number. 00305 * 00306 DO 120 I = 2, N - 1 00307 A( I ) = Y 00308 120 CONTINUE 00309 A( N ) = Z 00310 JC = N + 1 00311 DO 130 J = 2, N - 1 00312 A( JC+1 ) = WORK( J-1 ) 00313 IF( J.LT.N-1 ) 00314 $ A( JC+2 ) = WORK( N+J-1 ) 00315 A( JC+N-J ) = Y 00316 JC = JC + N - J + 1 00317 130 CONTINUE 00318 END IF 00319 * 00320 * Fill in the zeros using Givens rotations 00321 * 00322 IF( UPPER ) THEN 00323 JC = 1 00324 DO 150 J = 1, N - 1 00325 JCNEXT = JC + J 00326 RA = A( JCNEXT+J-1 ) 00327 RB = TWO 00328 CALL DROTG( RA, RB, C, S ) 00329 * 00330 * Multiply by [ c s; -s c] on the left. 00331 * 00332 IF( N.GT.J+1 ) THEN 00333 JX = JCNEXT + J 00334 DO 140 I = J + 2, N 00335 STEMP = C*A( JX+J ) + S*A( JX+J+1 ) 00336 A( JX+J+1 ) = -S*A( JX+J ) + C*A( JX+J+1 ) 00337 A( JX+J ) = STEMP 00338 JX = JX + I 00339 140 CONTINUE 00340 END IF 00341 * 00342 * Multiply by [-c -s; s -c] on the right. 00343 * 00344 IF( J.GT.1 ) 00345 $ CALL DROT( J-1, A( JCNEXT ), 1, A( JC ), 1, -C, -S ) 00346 * 00347 * Negate A(J,J+1). 00348 * 00349 A( JCNEXT+J-1 ) = -A( JCNEXT+J-1 ) 00350 JC = JCNEXT 00351 150 CONTINUE 00352 ELSE 00353 JC = 1 00354 DO 170 J = 1, N - 1 00355 JCNEXT = JC + N - J + 1 00356 RA = A( JC+1 ) 00357 RB = TWO 00358 CALL DROTG( RA, RB, C, S ) 00359 * 00360 * Multiply by [ c -s; s c] on the right. 00361 * 00362 IF( N.GT.J+1 ) 00363 $ CALL DROT( N-J-1, A( JCNEXT+1 ), 1, A( JC+2 ), 1, C, 00364 $ -S ) 00365 * 00366 * Multiply by [-c s; -s -c] on the left. 00367 * 00368 IF( J.GT.1 ) THEN 00369 JX = 1 00370 DO 160 I = 1, J - 1 00371 STEMP = -C*A( JX+J-I ) + S*A( JX+J-I+1 ) 00372 A( JX+J-I+1 ) = -S*A( JX+J-I ) - C*A( JX+J-I+1 ) 00373 A( JX+J-I ) = STEMP 00374 JX = JX + N - I + 1 00375 160 CONTINUE 00376 END IF 00377 * 00378 * Negate A(J+1,J). 00379 * 00380 A( JC+1 ) = -A( JC+1 ) 00381 JC = JCNEXT 00382 170 CONTINUE 00383 END IF 00384 * 00385 * IMAT > 10: Pathological test cases. These triangular matrices 00386 * are badly scaled or badly conditioned, so when used in solving a 00387 * triangular system they may cause overflow in the solution vector. 00388 * 00389 ELSE IF( IMAT.EQ.11 ) THEN 00390 * 00391 * Type 11: Generate a triangular matrix with elements between 00392 * -1 and 1. Give the diagonal norm 2 to make it well-conditioned. 00393 * Make the right hand side large so that it requires scaling. 00394 * 00395 IF( UPPER ) THEN 00396 JC = 1 00397 DO 180 J = 1, N 00398 CALL DLARNV( 2, ISEED, J, A( JC ) ) 00399 A( JC+J-1 ) = SIGN( TWO, A( JC+J-1 ) ) 00400 JC = JC + J 00401 180 CONTINUE 00402 ELSE 00403 JC = 1 00404 DO 190 J = 1, N 00405 CALL DLARNV( 2, ISEED, N-J+1, A( JC ) ) 00406 A( JC ) = SIGN( TWO, A( JC ) ) 00407 JC = JC + N - J + 1 00408 190 CONTINUE 00409 END IF 00410 * 00411 * Set the right hand side so that the largest value is BIGNUM. 00412 * 00413 CALL DLARNV( 2, ISEED, N, B ) 00414 IY = IDAMAX( N, B, 1 ) 00415 BNORM = ABS( B( IY ) ) 00416 BSCAL = BIGNUM / MAX( ONE, BNORM ) 00417 CALL DSCAL( N, BSCAL, B, 1 ) 00418 * 00419 ELSE IF( IMAT.EQ.12 ) THEN 00420 * 00421 * Type 12: Make the first diagonal element in the solve small to 00422 * cause immediate overflow when dividing by T(j,j). 00423 * In type 12, the offdiagonal elements are small (CNORM(j) < 1). 00424 * 00425 CALL DLARNV( 2, ISEED, N, B ) 00426 TSCAL = ONE / MAX( ONE, DBLE( N-1 ) ) 00427 IF( UPPER ) THEN 00428 JC = 1 00429 DO 200 J = 1, N 00430 CALL DLARNV( 2, ISEED, J-1, A( JC ) ) 00431 CALL DSCAL( J-1, TSCAL, A( JC ), 1 ) 00432 A( JC+J-1 ) = SIGN( ONE, DLARND( 2, ISEED ) ) 00433 JC = JC + J 00434 200 CONTINUE 00435 A( N*( N+1 ) / 2 ) = SMLNUM 00436 ELSE 00437 JC = 1 00438 DO 210 J = 1, N 00439 CALL DLARNV( 2, ISEED, N-J, A( JC+1 ) ) 00440 CALL DSCAL( N-J, TSCAL, A( JC+1 ), 1 ) 00441 A( JC ) = SIGN( ONE, DLARND( 2, ISEED ) ) 00442 JC = JC + N - J + 1 00443 210 CONTINUE 00444 A( 1 ) = SMLNUM 00445 END IF 00446 * 00447 ELSE IF( IMAT.EQ.13 ) THEN 00448 * 00449 * Type 13: Make the first diagonal element in the solve small to 00450 * cause immediate overflow when dividing by T(j,j). 00451 * In type 13, the offdiagonal elements are O(1) (CNORM(j) > 1). 00452 * 00453 CALL DLARNV( 2, ISEED, N, B ) 00454 IF( UPPER ) THEN 00455 JC = 1 00456 DO 220 J = 1, N 00457 CALL DLARNV( 2, ISEED, J-1, A( JC ) ) 00458 A( JC+J-1 ) = SIGN( ONE, DLARND( 2, ISEED ) ) 00459 JC = JC + J 00460 220 CONTINUE 00461 A( N*( N+1 ) / 2 ) = SMLNUM 00462 ELSE 00463 JC = 1 00464 DO 230 J = 1, N 00465 CALL DLARNV( 2, ISEED, N-J, A( JC+1 ) ) 00466 A( JC ) = SIGN( ONE, DLARND( 2, ISEED ) ) 00467 JC = JC + N - J + 1 00468 230 CONTINUE 00469 A( 1 ) = SMLNUM 00470 END IF 00471 * 00472 ELSE IF( IMAT.EQ.14 ) THEN 00473 * 00474 * Type 14: T is diagonal with small numbers on the diagonal to 00475 * make the growth factor underflow, but a small right hand side 00476 * chosen so that the solution does not overflow. 00477 * 00478 IF( UPPER ) THEN 00479 JCOUNT = 1 00480 JC = ( N-1 )*N / 2 + 1 00481 DO 250 J = N, 1, -1 00482 DO 240 I = 1, J - 1 00483 A( JC+I-1 ) = ZERO 00484 240 CONTINUE 00485 IF( JCOUNT.LE.2 ) THEN 00486 A( JC+J-1 ) = SMLNUM 00487 ELSE 00488 A( JC+J-1 ) = ONE 00489 END IF 00490 JCOUNT = JCOUNT + 1 00491 IF( JCOUNT.GT.4 ) 00492 $ JCOUNT = 1 00493 JC = JC - J + 1 00494 250 CONTINUE 00495 ELSE 00496 JCOUNT = 1 00497 JC = 1 00498 DO 270 J = 1, N 00499 DO 260 I = J + 1, N 00500 A( JC+I-J ) = ZERO 00501 260 CONTINUE 00502 IF( JCOUNT.LE.2 ) THEN 00503 A( JC ) = SMLNUM 00504 ELSE 00505 A( JC ) = ONE 00506 END IF 00507 JCOUNT = JCOUNT + 1 00508 IF( JCOUNT.GT.4 ) 00509 $ JCOUNT = 1 00510 JC = JC + N - J + 1 00511 270 CONTINUE 00512 END IF 00513 * 00514 * Set the right hand side alternately zero and small. 00515 * 00516 IF( UPPER ) THEN 00517 B( 1 ) = ZERO 00518 DO 280 I = N, 2, -2 00519 B( I ) = ZERO 00520 B( I-1 ) = SMLNUM 00521 280 CONTINUE 00522 ELSE 00523 B( N ) = ZERO 00524 DO 290 I = 1, N - 1, 2 00525 B( I ) = ZERO 00526 B( I+1 ) = SMLNUM 00527 290 CONTINUE 00528 END IF 00529 * 00530 ELSE IF( IMAT.EQ.15 ) THEN 00531 * 00532 * Type 15: Make the diagonal elements small to cause gradual 00533 * overflow when dividing by T(j,j). To control the amount of 00534 * scaling needed, the matrix is bidiagonal. 00535 * 00536 TEXP = ONE / MAX( ONE, DBLE( N-1 ) ) 00537 TSCAL = SMLNUM**TEXP 00538 CALL DLARNV( 2, ISEED, N, B ) 00539 IF( UPPER ) THEN 00540 JC = 1 00541 DO 310 J = 1, N 00542 DO 300 I = 1, J - 2 00543 A( JC+I-1 ) = ZERO 00544 300 CONTINUE 00545 IF( J.GT.1 ) 00546 $ A( JC+J-2 ) = -ONE 00547 A( JC+J-1 ) = TSCAL 00548 JC = JC + J 00549 310 CONTINUE 00550 B( N ) = ONE 00551 ELSE 00552 JC = 1 00553 DO 330 J = 1, N 00554 DO 320 I = J + 2, N 00555 A( JC+I-J ) = ZERO 00556 320 CONTINUE 00557 IF( J.LT.N ) 00558 $ A( JC+1 ) = -ONE 00559 A( JC ) = TSCAL 00560 JC = JC + N - J + 1 00561 330 CONTINUE 00562 B( 1 ) = ONE 00563 END IF 00564 * 00565 ELSE IF( IMAT.EQ.16 ) THEN 00566 * 00567 * Type 16: One zero diagonal element. 00568 * 00569 IY = N / 2 + 1 00570 IF( UPPER ) THEN 00571 JC = 1 00572 DO 340 J = 1, N 00573 CALL DLARNV( 2, ISEED, J, A( JC ) ) 00574 IF( J.NE.IY ) THEN 00575 A( JC+J-1 ) = SIGN( TWO, A( JC+J-1 ) ) 00576 ELSE 00577 A( JC+J-1 ) = ZERO 00578 END IF 00579 JC = JC + J 00580 340 CONTINUE 00581 ELSE 00582 JC = 1 00583 DO 350 J = 1, N 00584 CALL DLARNV( 2, ISEED, N-J+1, A( JC ) ) 00585 IF( J.NE.IY ) THEN 00586 A( JC ) = SIGN( TWO, A( JC ) ) 00587 ELSE 00588 A( JC ) = ZERO 00589 END IF 00590 JC = JC + N - J + 1 00591 350 CONTINUE 00592 END IF 00593 CALL DLARNV( 2, ISEED, N, B ) 00594 CALL DSCAL( N, TWO, B, 1 ) 00595 * 00596 ELSE IF( IMAT.EQ.17 ) THEN 00597 * 00598 * Type 17: Make the offdiagonal elements large to cause overflow 00599 * when adding a column of T. In the non-transposed case, the 00600 * matrix is constructed to cause overflow when adding a column in 00601 * every other step. 00602 * 00603 TSCAL = UNFL / ULP 00604 TSCAL = ( ONE-ULP ) / TSCAL 00605 DO 360 J = 1, N*( N+1 ) / 2 00606 A( J ) = ZERO 00607 360 CONTINUE 00608 TEXP = ONE 00609 IF( UPPER ) THEN 00610 JC = ( N-1 )*N / 2 + 1 00611 DO 370 J = N, 2, -2 00612 A( JC ) = -TSCAL / DBLE( N+1 ) 00613 A( JC+J-1 ) = ONE 00614 B( J ) = TEXP*( ONE-ULP ) 00615 JC = JC - J + 1 00616 A( JC ) = -( TSCAL / DBLE( N+1 ) ) / DBLE( N+2 ) 00617 A( JC+J-2 ) = ONE 00618 B( J-1 ) = TEXP*DBLE( N*N+N-1 ) 00619 TEXP = TEXP*TWO 00620 JC = JC - J + 2 00621 370 CONTINUE 00622 B( 1 ) = ( DBLE( N+1 ) / DBLE( N+2 ) )*TSCAL 00623 ELSE 00624 JC = 1 00625 DO 380 J = 1, N - 1, 2 00626 A( JC+N-J ) = -TSCAL / DBLE( N+1 ) 00627 A( JC ) = ONE 00628 B( J ) = TEXP*( ONE-ULP ) 00629 JC = JC + N - J + 1 00630 A( JC+N-J-1 ) = -( TSCAL / DBLE( N+1 ) ) / DBLE( N+2 ) 00631 A( JC ) = ONE 00632 B( J+1 ) = TEXP*DBLE( N*N+N-1 ) 00633 TEXP = TEXP*TWO 00634 JC = JC + N - J 00635 380 CONTINUE 00636 B( N ) = ( DBLE( N+1 ) / DBLE( N+2 ) )*TSCAL 00637 END IF 00638 * 00639 ELSE IF( IMAT.EQ.18 ) THEN 00640 * 00641 * Type 18: Generate a unit triangular matrix with elements 00642 * between -1 and 1, and make the right hand side large so that it 00643 * requires scaling. 00644 * 00645 IF( UPPER ) THEN 00646 JC = 1 00647 DO 390 J = 1, N 00648 CALL DLARNV( 2, ISEED, J-1, A( JC ) ) 00649 A( JC+J-1 ) = ZERO 00650 JC = JC + J 00651 390 CONTINUE 00652 ELSE 00653 JC = 1 00654 DO 400 J = 1, N 00655 IF( J.LT.N ) 00656 $ CALL DLARNV( 2, ISEED, N-J, A( JC+1 ) ) 00657 A( JC ) = ZERO 00658 JC = JC + N - J + 1 00659 400 CONTINUE 00660 END IF 00661 * 00662 * Set the right hand side so that the largest value is BIGNUM. 00663 * 00664 CALL DLARNV( 2, ISEED, N, B ) 00665 IY = IDAMAX( N, B, 1 ) 00666 BNORM = ABS( B( IY ) ) 00667 BSCAL = BIGNUM / MAX( ONE, BNORM ) 00668 CALL DSCAL( N, BSCAL, B, 1 ) 00669 * 00670 ELSE IF( IMAT.EQ.19 ) THEN 00671 * 00672 * Type 19: Generate a triangular matrix with elements between 00673 * BIGNUM/(n-1) and BIGNUM so that at least one of the column 00674 * norms will exceed BIGNUM. 00675 * 00676 TLEFT = BIGNUM / MAX( ONE, DBLE( N-1 ) ) 00677 TSCAL = BIGNUM*( DBLE( N-1 ) / MAX( ONE, DBLE( N ) ) ) 00678 IF( UPPER ) THEN 00679 JC = 1 00680 DO 420 J = 1, N 00681 CALL DLARNV( 2, ISEED, J, A( JC ) ) 00682 DO 410 I = 1, J 00683 A( JC+I-1 ) = SIGN( TLEFT, A( JC+I-1 ) ) + 00684 $ TSCAL*A( JC+I-1 ) 00685 410 CONTINUE 00686 JC = JC + J 00687 420 CONTINUE 00688 ELSE 00689 JC = 1 00690 DO 440 J = 1, N 00691 CALL DLARNV( 2, ISEED, N-J+1, A( JC ) ) 00692 DO 430 I = J, N 00693 A( JC+I-J ) = SIGN( TLEFT, A( JC+I-J ) ) + 00694 $ TSCAL*A( JC+I-J ) 00695 430 CONTINUE 00696 JC = JC + N - J + 1 00697 440 CONTINUE 00698 END IF 00699 CALL DLARNV( 2, ISEED, N, B ) 00700 CALL DSCAL( N, TWO, B, 1 ) 00701 END IF 00702 * 00703 * Flip the matrix across its counter-diagonal if the transpose will 00704 * be used. 00705 * 00706 IF( .NOT.LSAME( TRANS, 'N' ) ) THEN 00707 IF( UPPER ) THEN 00708 JJ = 1 00709 JR = N*( N+1 ) / 2 00710 DO 460 J = 1, N / 2 00711 JL = JJ 00712 DO 450 I = J, N - J 00713 T = A( JR-I+J ) 00714 A( JR-I+J ) = A( JL ) 00715 A( JL ) = T 00716 JL = JL + I 00717 450 CONTINUE 00718 JJ = JJ + J + 1 00719 JR = JR - ( N-J+1 ) 00720 460 CONTINUE 00721 ELSE 00722 JL = 1 00723 JJ = N*( N+1 ) / 2 00724 DO 480 J = 1, N / 2 00725 JR = JJ 00726 DO 470 I = J, N - J 00727 T = A( JL+I-J ) 00728 A( JL+I-J ) = A( JR ) 00729 A( JR ) = T 00730 JR = JR - I 00731 470 CONTINUE 00732 JL = JL + N - J + 1 00733 JJ = JJ - J - 1 00734 480 CONTINUE 00735 END IF 00736 END IF 00737 * 00738 RETURN 00739 * 00740 * End of DLATTP 00741 * 00742 END
1.7.3 
