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