LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlattp.f
Go to the documentation of this file.
00001 *> \brief \b DLATTP
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 DLATTP( 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 *       DOUBLE PRECISION   A( * ), B( * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DLATTP 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 *>          DLATMS).  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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
00097 *>          The right hand side vector, if IMAT > 10.
00098 *> \endverbatim
00099 *>
00100 *> \param[out] WORK
00101 *> \verbatim
00102 *>          WORK is DOUBLE PRECISION 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 double_lin
00123 *
00124 *  =====================================================================
00125       SUBROUTINE DLATTP( 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       DOUBLE PRECISION   A( * ), B( * ), WORK( * )
00140 *     ..
00141 *
00142 *  =====================================================================
00143 *
00144 *     .. Parameters ..
00145       DOUBLE PRECISION   ONE, TWO, ZERO
00146       PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0, ZERO = 0.0D+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       DOUBLE PRECISION   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            IDAMAX
00162       DOUBLE PRECISION   DLAMCH, DLARND
00163       EXTERNAL           LSAME, IDAMAX, DLAMCH, DLARND
00164 *     ..
00165 *     .. External Subroutines ..
00166       EXTERNAL           DLABAD, DLARNV, DLATB4, DLATMS, DROT, DROTG,
00167      $                   DSCAL
00168 *     ..
00169 *     .. Intrinsic Functions ..
00170       INTRINSIC          ABS, DBLE, MAX, SIGN, SQRT
00171 *     ..
00172 *     .. Executable Statements ..
00173 *
00174       PATH( 1: 1 ) = 'Double precision'
00175       PATH( 2: 3 ) = 'TP'
00176       UNFL = DLAMCH( 'Safe minimum' )
00177       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00178       SMLNUM = UNFL
00179       BIGNUM = ( ONE-ULP ) / SMLNUM
00180       CALL DLABAD( 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 DLATB4 to set parameters for SLATMS.
00194 *
00195       UPPER = LSAME( UPLO, 'U' )
00196       IF( UPPER ) THEN
00197          CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00198      $                CNDNUM, DIST )
00199          PACKIT = 'C'
00200       ELSE
00201          CALL DLATB4( 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 DLATMS( 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.25D0
00324          SFAC = 0.5D0
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 = DLARND( 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 / DBLE( 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 DROTG( 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 DROT( 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 DROTG( RA, RB, C, S )
00429 *
00430 *              Multiply by [ c -s;  s  c] on the right.
00431 *
00432                IF( N.GT.J+1 )
00433      $            CALL DROT( 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 DLARNV( 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 DLARNV( 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 DLARNV( 2, ISEED, N, B )
00484          IY = IDAMAX( N, B, 1 )
00485          BNORM = ABS( B( IY ) )
00486          BSCAL = BIGNUM / MAX( ONE, BNORM )
00487          CALL DSCAL( 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 DLARNV( 2, ISEED, N, B )
00496          TSCAL = ONE / MAX( ONE, DBLE( N-1 ) )
00497          IF( UPPER ) THEN
00498             JC = 1
00499             DO 200 J = 1, N
00500                CALL DLARNV( 2, ISEED, J-1, A( JC ) )
00501                CALL DSCAL( J-1, TSCAL, A( JC ), 1 )
00502                A( JC+J-1 ) = SIGN( ONE, DLARND( 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 DLARNV( 2, ISEED, N-J, A( JC+1 ) )
00510                CALL DSCAL( N-J, TSCAL, A( JC+1 ), 1 )
00511                A( JC ) = SIGN( ONE, DLARND( 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 DLARNV( 2, ISEED, N, B )
00524          IF( UPPER ) THEN
00525             JC = 1
00526             DO 220 J = 1, N
00527                CALL DLARNV( 2, ISEED, J-1, A( JC ) )
00528                A( JC+J-1 ) = SIGN( ONE, DLARND( 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 DLARNV( 2, ISEED, N-J, A( JC+1 ) )
00536                A( JC ) = SIGN( ONE, DLARND( 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, DBLE( N-1 ) )
00607          TSCAL = SMLNUM**TEXP
00608          CALL DLARNV( 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 DLARNV( 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 DLARNV( 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 DLARNV( 2, ISEED, N, B )
00664          CALL DSCAL( 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 / DBLE( N+1 )
00683                A( JC+J-1 ) = ONE
00684                B( J ) = TEXP*( ONE-ULP )
00685                JC = JC - J + 1
00686                A( JC ) = -( TSCAL / DBLE( N+1 ) ) / DBLE( N+2 )
00687                A( JC+J-2 ) = ONE
00688                B( J-1 ) = TEXP*DBLE( N*N+N-1 )
00689                TEXP = TEXP*TWO
00690                JC = JC - J + 2
00691   370       CONTINUE
00692             B( 1 ) = ( DBLE( N+1 ) / DBLE( N+2 ) )*TSCAL
00693          ELSE
00694             JC = 1
00695             DO 380 J = 1, N - 1, 2
00696                A( JC+N-J ) = -TSCAL / DBLE( 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 / DBLE( N+1 ) ) / DBLE( N+2 )
00701                A( JC ) = ONE
00702                B( J+1 ) = TEXP*DBLE( N*N+N-1 )
00703                TEXP = TEXP*TWO
00704                JC = JC + N - J
00705   380       CONTINUE
00706             B( N ) = ( DBLE( N+1 ) / DBLE( 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 DLARNV( 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 DLARNV( 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 DLARNV( 2, ISEED, N, B )
00735          IY = IDAMAX( N, B, 1 )
00736          BNORM = ABS( B( IY ) )
00737          BSCAL = BIGNUM / MAX( ONE, BNORM )
00738          CALL DSCAL( 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, DBLE( N-1 ) )
00747          TSCAL = BIGNUM*( DBLE( N-1 ) / MAX( ONE, DBLE( N ) ) )
00748          IF( UPPER ) THEN
00749             JC = 1
00750             DO 420 J = 1, N
00751                CALL DLARNV( 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 DLARNV( 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 DLARNV( 2, ISEED, N, B )
00770          CALL DSCAL( 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 DLATTP
00811 *
00812       END
 All Files Functions