LAPACK 3.3.0

slattp.f

Go to the documentation of this file.
00001       SUBROUTINE SLATTP( 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       REAL               A( * ), B( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SLATTP 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 *          SLATMS).  Modified on exit.
00051 *
00052 *  N       (input) INTEGER
00053 *          The order of the matrix to be generated.
00054 *
00055 *  A       (output) REAL 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) REAL array, dimension (N)
00064 *          The right hand side vector, if IMAT > 10.
00065 *
00066 *  WORK    (workspace) REAL 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       REAL               ONE, TWO, ZERO
00076       PARAMETER          ( ONE = 1.0E+0, TWO = 2.0E+0, ZERO = 0.0E+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       REAL               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            ISAMAX
00092       REAL               SLAMCH, SLARND
00093       EXTERNAL           LSAME, ISAMAX, SLAMCH, SLARND
00094 *     ..
00095 *     .. External Subroutines ..
00096       EXTERNAL           SLABAD, SLARNV, SLATB4, SLATMS, SROT, SROTG,
00097      $                   SSCAL
00098 *     ..
00099 *     .. Intrinsic Functions ..
00100       INTRINSIC          ABS, MAX, REAL, SIGN, SQRT
00101 *     ..
00102 *     .. Executable Statements ..
00103 *
00104       PATH( 1: 1 ) = 'Single precision'
00105       PATH( 2: 3 ) = 'TP'
00106       UNFL = SLAMCH( 'Safe minimum' )
00107       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00108       SMLNUM = UNFL
00109       BIGNUM = ( ONE-ULP ) / SMLNUM
00110       CALL SLABAD( 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 SLATB4 to set parameters for SLATMS.
00124 *
00125       UPPER = LSAME( UPLO, 'U' )
00126       IF( UPPER ) THEN
00127          CALL SLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00128      $                CNDNUM, DIST )
00129          PACKIT = 'C'
00130       ELSE
00131          CALL SLATB4( 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 SLATMS( 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.25
00254          SFAC = 0.5
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 = SLARND( 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 / REAL( 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 SROTG( 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 SROT( 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 SROTG( RA, RB, C, S )
00359 *
00360 *              Multiply by [ c -s;  s  c] on the right.
00361 *
00362                IF( N.GT.J+1 )
00363      $            CALL SROT( 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 SLARNV( 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 SLARNV( 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 SLARNV( 2, ISEED, N, B )
00414          IY = ISAMAX( N, B, 1 )
00415          BNORM = ABS( B( IY ) )
00416          BSCAL = BIGNUM / MAX( ONE, BNORM )
00417          CALL SSCAL( 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 SLARNV( 2, ISEED, N, B )
00426          TSCAL = ONE / MAX( ONE, REAL( N-1 ) )
00427          IF( UPPER ) THEN
00428             JC = 1
00429             DO 200 J = 1, N
00430                CALL SLARNV( 2, ISEED, J-1, A( JC ) )
00431                CALL SSCAL( J-1, TSCAL, A( JC ), 1 )
00432                A( JC+J-1 ) = SIGN( ONE, SLARND( 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 SLARNV( 2, ISEED, N-J, A( JC+1 ) )
00440                CALL SSCAL( N-J, TSCAL, A( JC+1 ), 1 )
00441                A( JC ) = SIGN( ONE, SLARND( 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 SLARNV( 2, ISEED, N, B )
00454          IF( UPPER ) THEN
00455             JC = 1
00456             DO 220 J = 1, N
00457                CALL SLARNV( 2, ISEED, J-1, A( JC ) )
00458                A( JC+J-1 ) = SIGN( ONE, SLARND( 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 SLARNV( 2, ISEED, N-J, A( JC+1 ) )
00466                A( JC ) = SIGN( ONE, SLARND( 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, REAL( N-1 ) )
00537          TSCAL = SMLNUM**TEXP
00538          CALL SLARNV( 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 SLARNV( 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 SLARNV( 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 SLARNV( 2, ISEED, N, B )
00594          CALL SSCAL( 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 / REAL( N+1 )
00613                A( JC+J-1 ) = ONE
00614                B( J ) = TEXP*( ONE-ULP )
00615                JC = JC - J + 1
00616                A( JC ) = -( TSCAL / REAL( N+1 ) ) / REAL( N+2 )
00617                A( JC+J-2 ) = ONE
00618                B( J-1 ) = TEXP*REAL( N*N+N-1 )
00619                TEXP = TEXP*TWO
00620                JC = JC - J + 2
00621   370       CONTINUE
00622             B( 1 ) = ( REAL( N+1 ) / REAL( N+2 ) )*TSCAL
00623          ELSE
00624             JC = 1
00625             DO 380 J = 1, N - 1, 2
00626                A( JC+N-J ) = -TSCAL / REAL( 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 / REAL( N+1 ) ) / REAL( N+2 )
00631                A( JC ) = ONE
00632                B( J+1 ) = TEXP*REAL( N*N+N-1 )
00633                TEXP = TEXP*TWO
00634                JC = JC + N - J
00635   380       CONTINUE
00636             B( N ) = ( REAL( N+1 ) / REAL( 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 SLARNV( 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 SLARNV( 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 SLARNV( 2, ISEED, N, B )
00665          IY = ISAMAX( N, B, 1 )
00666          BNORM = ABS( B( IY ) )
00667          BSCAL = BIGNUM / MAX( ONE, BNORM )
00668          CALL SSCAL( 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, REAL( N-1 ) )
00677          TSCAL = BIGNUM*( REAL( N-1 ) / MAX( ONE, REAL( N ) ) )
00678          IF( UPPER ) THEN
00679             JC = 1
00680             DO 420 J = 1, N
00681                CALL SLARNV( 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 SLARNV( 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 SLARNV( 2, ISEED, N, B )
00700          CALL SSCAL( 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 SLATTP
00741 *
00742       END
 All Files Functions