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