LAPACK: TESTING/EIG/dchksb.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DCHKSB( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, ISEED,
00002      $                   THRESH, NOUNIT, A, LDA, SD, SE, U, LDU, WORK,
00003      $                   LWORK, RESULT, INFO )
00004 *
00005 *  -- LAPACK test routine (version 3.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
00011      $                   NWDTHS
00012       DOUBLE PRECISION   THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            ISEED( 4 ), KK( * ), NN( * )
00017       DOUBLE PRECISION   A( LDA, * ), RESULT( * ), SD( * ), SE( * ),
00018      $                   U( LDU, * ), WORK( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  DCHKSB tests the reduction of a symmetric band matrix to tridiagonal
00025 *  form, used with the symmetric eigenvalue problem.
00026 *
00027 *  DSBTRD factors a symmetric band matrix A as  U S U' , where ' means
00028 *  transpose, S is symmetric tridiagonal, and U is orthogonal.
00029 *  DSBTRD can use either just the lower or just the upper triangle
00030 *  of A; DCHKSB checks both cases.
00031 *
00032 *  When DCHKSB is called, a number of matrix "sizes" ("n's"), a number
00033 *  of bandwidths ("k's"), and a number of matrix "types" are
00034 *  specified.  For each size ("n"), each bandwidth ("k") less than or
00035 *  equal to "n", and each type of matrix, one matrix will be generated
00036 *  and used to test the symmetric banded reduction routine.  For each
00037 *  matrix, a number of tests will be performed:
00038 *
00039 *  (1)     | A - V S V' | / ( |A| n ulp )  computed by DSBTRD with
00040 *                                          UPLO='U'
00041 *
00042 *  (2)     | I - UU' | / ( n ulp )
00043 *
00044 *  (3)     | A - V S V' | / ( |A| n ulp )  computed by DSBTRD with
00045 *                                          UPLO='L'
00046 *
00047 *  (4)     | I - UU' | / ( n ulp )
00048 *
00049 *  The "sizes" are specified by an array NN(1:NSIZES); the value of
00050 *  each element NN(j) specifies one size.
00051 *  The "types" are specified by a logical array DOTYPE( 1:NTYPES );
00052 *  if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
00053 *  Currently, the list of possible types is:
00054 *
00055 *  (1)  The zero matrix.
00056 *  (2)  The identity matrix.
00057 *
00058 *  (3)  A diagonal matrix with evenly spaced entries
00059 *       1, ..., ULP  and random signs.
00060 *       (ULP = (first number larger than 1) - 1 )
00061 *  (4)  A diagonal matrix with geometrically spaced entries
00062 *       1, ..., ULP  and random signs.
00063 *  (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
00064 *       and random signs.
00065 *
00066 *  (6)  Same as (4), but multiplied by SQRT( overflow threshold )
00067 *  (7)  Same as (4), but multiplied by SQRT( underflow threshold )
00068 *
00069 *  (8)  A matrix of the form  U' D U, where U is orthogonal and
00070 *       D has evenly spaced entries 1, ..., ULP with random signs
00071 *       on the diagonal.
00072 *
00073 *  (9)  A matrix of the form  U' D U, where U is orthogonal and
00074 *       D has geometrically spaced entries 1, ..., ULP with random
00075 *       signs on the diagonal.
00076 *
00077 *  (10) A matrix of the form  U' D U, where U is orthogonal and
00078 *       D has "clustered" entries 1, ULP,..., ULP with random
00079 *       signs on the diagonal.
00080 *
00081 *  (11) Same as (8), but multiplied by SQRT( overflow threshold )
00082 *  (12) Same as (8), but multiplied by SQRT( underflow threshold )
00083 *
00084 *  (13) Symmetric matrix with random entries chosen from (-1,1).
00085 *  (14) Same as (13), but multiplied by SQRT( overflow threshold )
00086 *  (15) Same as (13), but multiplied by SQRT( underflow threshold )
00087 *
00088 *  Arguments
00089 *  =========
00090 *
00091 *  NSIZES  (input) INTEGER
00092 *          The number of sizes of matrices to use.  If it is zero,
00093 *          DCHKSB does nothing.  It must be at least zero.
00094 *
00095 *  NN      (input) INTEGER array, dimension (NSIZES)
00096 *          An array containing the sizes to be used for the matrices.
00097 *          Zero values will be skipped.  The values must be at least
00098 *          zero.
00099 *
00100 *  NWDTHS  (input) INTEGER
00101 *          The number of bandwidths to use.  If it is zero,
00102 *          DCHKSB does nothing.  It must be at least zero.
00103 *
00104 *  KK      (input) INTEGER array, dimension (NWDTHS)
00105 *          An array containing the bandwidths to be used for the band
00106 *          matrices.  The values must be at least zero.
00107 *
00108 *  NTYPES  (input) INTEGER
00109 *          The number of elements in DOTYPE.   If it is zero, DCHKSB
00110 *          does nothing.  It must be at least zero.  If it is MAXTYP+1
00111 *          and NSIZES is 1, then an additional type, MAXTYP+1 is
00112 *          defined, which is to use whatever matrix is in A.  This
00113 *          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00114 *          DOTYPE(MAXTYP+1) is .TRUE. .
00115 *
00116 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00117 *          If DOTYPE(j) is .TRUE., then for each size in NN a
00118 *          matrix of that size and of type j will be generated.
00119 *          If NTYPES is smaller than the maximum number of types
00120 *          defined (PARAMETER MAXTYP), then types NTYPES+1 through
00121 *          MAXTYP will not be generated.  If NTYPES is larger
00122 *          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
00123 *          will be ignored.
00124 *
00125 *  ISEED   (input/output) INTEGER array, dimension (4)
00126 *          On entry ISEED specifies the seed of the random number
00127 *          generator. The array elements should be between 0 and 4095;
00128 *          if not they will be reduced mod 4096.  Also, ISEED(4) must
00129 *          be odd.  The random number generator uses a linear
00130 *          congruential sequence limited to small integers, and so
00131 *          should produce machine independent random numbers. The
00132 *          values of ISEED are changed on exit, and can be used in the
00133 *          next call to DCHKSB to continue the same random number
00134 *          sequence.
00135 *
00136 *  THRESH  (input) DOUBLE PRECISION
00137 *          A test will count as "failed" if the "error", computed as
00138 *          described above, exceeds THRESH.  Note that the error
00139 *          is scaled to be O(1), so THRESH should be a reasonably
00140 *          small multiple of 1, e.g., 10 or 100.  In particular,
00141 *          it should not depend on the precision (single vs. double)
00142 *          or the size of the matrix.  It must be at least zero.
00143 *
00144 *  NOUNIT  (input) INTEGER
00145 *          The FORTRAN unit number for printing out error messages
00146 *          (e.g., if a routine returns IINFO not equal to 0.)
00147 *
00148 *  A       (input/workspace) DOUBLE PRECISION array, dimension
00149 *                            (LDA, max(NN))
00150 *          Used to hold the matrix whose eigenvalues are to be
00151 *          computed.
00152 *
00153 *  LDA     (input) INTEGER
00154 *          The leading dimension of A.  It must be at least 2 (not 1!)
00155 *          and at least max( KK )+1.
00156 *
00157 *  SD      (workspace) DOUBLE PRECISION array, dimension (max(NN))
00158 *          Used to hold the diagonal of the tridiagonal matrix computed
00159 *          by DSBTRD.
00160 *
00161 *  SE      (workspace) DOUBLE PRECISION array, dimension (max(NN))
00162 *          Used to hold the off-diagonal of the tridiagonal matrix
00163 *          computed by DSBTRD.
00164 *
00165 *  U       (workspace) DOUBLE PRECISION array, dimension (LDU, max(NN))
00166 *          Used to hold the orthogonal matrix computed by DSBTRD.
00167 *
00168 *  LDU     (input) INTEGER
00169 *          The leading dimension of U.  It must be at least 1
00170 *          and at least max( NN ).
00171 *
00172 *  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
00173 *
00174 *  LWORK   (input) INTEGER
00175 *          The number of entries in WORK.  This must be at least
00176 *          max( LDA+1, max(NN)+1 )*max(NN).
00177 *
00178 *  RESULT  (output) DOUBLE PRECISION array, dimension (4)
00179 *          The values computed by the tests described above.
00180 *          The values are currently limited to 1/ulp, to avoid
00181 *          overflow.
00182 *
00183 *  INFO    (output) INTEGER
00184 *          If 0, then everything ran OK.
00185 *
00186 *-----------------------------------------------------------------------
00187 *
00188 *       Some Local Variables and Parameters:
00189 *       ---- ----- --------- --- ----------
00190 *       ZERO, ONE       Real 0 and 1.
00191 *       MAXTYP          The number of types defined.
00192 *       NTEST           The number of tests performed, or which can
00193 *                       be performed so far, for the current matrix.
00194 *       NTESTT          The total number of tests performed so far.
00195 *       NMAX            Largest value in NN.
00196 *       NMATS           The number of matrices generated so far.
00197 *       NERRS           The number of tests which have exceeded THRESH
00198 *                       so far.
00199 *       COND, IMODE     Values to be passed to the matrix generators.
00200 *       ANORM           Norm of A; passed to matrix generators.
00201 *
00202 *       OVFL, UNFL      Overflow and underflow thresholds.
00203 *       ULP, ULPINV     Finest relative precision and its inverse.
00204 *       RTOVFL, RTUNFL  Square roots of the previous 2 values.
00205 *               The following four arrays decode JTYPE:
00206 *       KTYPE(j)        The general type (1-10) for type "j".
00207 *       KMODE(j)        The MODE value to be passed to the matrix
00208 *                       generator for type "j".
00209 *       KMAGN(j)        The order of magnitude ( O(1),
00210 *                       O(overflow^(1/2) ), O(underflow^(1/2) )
00211 *
00212 *  =====================================================================
00213 *
00214 *     .. Parameters ..
00215       DOUBLE PRECISION   ZERO, ONE, TWO, TEN
00216       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00217      $                   TEN = 10.0D0 )
00218       DOUBLE PRECISION   HALF
00219       PARAMETER          ( HALF = ONE / TWO )
00220       INTEGER            MAXTYP
00221       PARAMETER          ( MAXTYP = 15 )
00222 *     ..
00223 *     .. Local Scalars ..
00224       LOGICAL            BADNN, BADNNB
00225       INTEGER            I, IINFO, IMODE, ITYPE, J, JC, JCOL, JR, JSIZE,
00226      $                   JTYPE, JWIDTH, K, KMAX, MTYPES, N, NERRS,
00227      $                   NMATS, NMAX, NTEST, NTESTT
00228       DOUBLE PRECISION   ANINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
00229      $                   TEMP1, ULP, ULPINV, UNFL
00230 *     ..
00231 *     .. Local Arrays ..
00232       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ),
00233      $                   KMODE( MAXTYP ), KTYPE( MAXTYP )
00234 *     ..
00235 *     .. External Functions ..
00236       DOUBLE PRECISION   DLAMCH
00237       EXTERNAL           DLAMCH
00238 *     ..
00239 *     .. External Subroutines ..
00240       EXTERNAL           DLACPY, DLASET, DLASUM, DLATMR, DLATMS, DSBT21,
00241      $                   DSBTRD, XERBLA
00242 *     ..
00243 *     .. Intrinsic Functions ..
00244       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
00245 *     ..
00246 *     .. Data statements ..
00247       DATA               KTYPE / 1, 2, 5*4, 5*5, 3*8 /
00248       DATA               KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
00249      $                   2, 3 /
00250       DATA               KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
00251      $                   0, 0 /
00252 *     ..
00253 *     .. Executable Statements ..
00254 *
00255 *     Check for errors
00256 *
00257       NTESTT = 0
00258       INFO = 0
00259 *
00260 *     Important constants
00261 *
00262       BADNN = .FALSE.
00263       NMAX = 1
00264       DO 10 J = 1, NSIZES
00265          NMAX = MAX( NMAX, NN( J ) )
00266          IF( NN( J ).LT.0 )
00267      $      BADNN = .TRUE.
00268    10 CONTINUE
00269 *
00270       BADNNB = .FALSE.
00271       KMAX = 0
00272       DO 20 J = 1, NSIZES
00273          KMAX = MAX( KMAX, KK( J ) )
00274          IF( KK( J ).LT.0 )
00275      $      BADNNB = .TRUE.
00276    20 CONTINUE
00277       KMAX = MIN( NMAX-1, KMAX )
00278 *
00279 *     Check for errors
00280 *
00281       IF( NSIZES.LT.0 ) THEN
00282          INFO = -1
00283       ELSE IF( BADNN ) THEN
00284          INFO = -2
00285       ELSE IF( NWDTHS.LT.0 ) THEN
00286          INFO = -3
00287       ELSE IF( BADNNB ) THEN
00288          INFO = -4
00289       ELSE IF( NTYPES.LT.0 ) THEN
00290          INFO = -5
00291       ELSE IF( LDA.LT.KMAX+1 ) THEN
00292          INFO = -11
00293       ELSE IF( LDU.LT.NMAX ) THEN
00294          INFO = -15
00295       ELSE IF( ( MAX( LDA, NMAX )+1 )*NMAX.GT.LWORK ) THEN
00296          INFO = -17
00297       END IF
00298 *
00299       IF( INFO.NE.0 ) THEN
00300          CALL XERBLA( 'DCHKSB', -INFO )
00301          RETURN
00302       END IF
00303 *
00304 *     Quick return if possible
00305 *
00306       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 .OR. NWDTHS.EQ.0 )
00307      $   RETURN
00308 *
00309 *     More Important constants
00310 *
00311       UNFL = DLAMCH( 'Safe minimum' )
00312       OVFL = ONE / UNFL
00313       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00314       ULPINV = ONE / ULP
00315       RTUNFL = SQRT( UNFL )
00316       RTOVFL = SQRT( OVFL )
00317 *
00318 *     Loop over sizes, types
00319 *
00320       NERRS = 0
00321       NMATS = 0
00322 *
00323       DO 190 JSIZE = 1, NSIZES
00324          N = NN( JSIZE )
00325          ANINV = ONE / DBLE( MAX( 1, N ) )
00326 *
00327          DO 180 JWIDTH = 1, NWDTHS
00328             K = KK( JWIDTH )
00329             IF( K.GT.N )
00330      $         GO TO 180
00331             K = MAX( 0, MIN( N-1, K ) )
00332 *
00333             IF( NSIZES.NE.1 ) THEN
00334                MTYPES = MIN( MAXTYP, NTYPES )
00335             ELSE
00336                MTYPES = MIN( MAXTYP+1, NTYPES )
00337             END IF
00338 *
00339             DO 170 JTYPE = 1, MTYPES
00340                IF( .NOT.DOTYPE( JTYPE ) )
00341      $            GO TO 170
00342                NMATS = NMATS + 1
00343                NTEST = 0
00344 *
00345                DO 30 J = 1, 4
00346                   IOLDSD( J ) = ISEED( J )
00347    30          CONTINUE
00348 *
00349 *              Compute "A".
00350 *              Store as "Upper"; later, we will copy to other format.
00351 *
00352 *              Control parameters:
00353 *
00354 *                  KMAGN  KMODE        KTYPE
00355 *              =1  O(1)   clustered 1  zero
00356 *              =2  large  clustered 2  identity
00357 *              =3  small  exponential  (none)
00358 *              =4         arithmetic   diagonal, (w/ eigenvalues)
00359 *              =5         random log   symmetric, w/ eigenvalues
00360 *              =6         random       (none)
00361 *              =7                      random diagonal
00362 *              =8                      random symmetric
00363 *              =9                      positive definite
00364 *              =10                     diagonally dominant tridiagonal
00365 *
00366                IF( MTYPES.GT.MAXTYP )
00367      $            GO TO 100
00368 *
00369                ITYPE = KTYPE( JTYPE )
00370                IMODE = KMODE( JTYPE )
00371 *
00372 *              Compute norm
00373 *
00374                GO TO ( 40, 50, 60 )KMAGN( JTYPE )
00375 *
00376    40          CONTINUE
00377                ANORM = ONE
00378                GO TO 70
00379 *
00380    50          CONTINUE
00381                ANORM = ( RTOVFL*ULP )*ANINV
00382                GO TO 70
00383 *
00384    60          CONTINUE
00385                ANORM = RTUNFL*N*ULPINV
00386                GO TO 70
00387 *
00388    70          CONTINUE
00389 *
00390                CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
00391                IINFO = 0
00392                IF( JTYPE.LE.15 ) THEN
00393                   COND = ULPINV
00394                ELSE
00395                   COND = ULPINV*ANINV / TEN
00396                END IF
00397 *
00398 *              Special Matrices -- Identity & Jordan block
00399 *
00400 *                 Zero
00401 *
00402                IF( ITYPE.EQ.1 ) THEN
00403                   IINFO = 0
00404 *
00405                ELSE IF( ITYPE.EQ.2 ) THEN
00406 *
00407 *                 Identity
00408 *
00409                   DO 80 JCOL = 1, N
00410                      A( K+1, JCOL ) = ANORM
00411    80             CONTINUE
00412 *
00413                ELSE IF( ITYPE.EQ.4 ) THEN
00414 *
00415 *                 Diagonal Matrix, [Eigen]values Specified
00416 *
00417                   CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
00418      $                         ANORM, 0, 0, 'Q', A( K+1, 1 ), LDA,
00419      $                         WORK( N+1 ), IINFO )
00420 *
00421                ELSE IF( ITYPE.EQ.5 ) THEN
00422 *
00423 *                 Symmetric, eigenvalues specified
00424 *
00425                   CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
00426      $                         ANORM, K, K, 'Q', A, LDA, WORK( N+1 ),
00427      $                         IINFO )
00428 *
00429                ELSE IF( ITYPE.EQ.7 ) THEN
00430 *
00431 *                 Diagonal, random eigenvalues
00432 *
00433                   CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
00434      $                         'T', 'N', WORK( N+1 ), 1, ONE,
00435      $                         WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00436      $                         ZERO, ANORM, 'Q', A( K+1, 1 ), LDA,
00437      $                         IDUMMA, IINFO )
00438 *
00439                ELSE IF( ITYPE.EQ.8 ) THEN
00440 *
00441 *                 Symmetric, random eigenvalues
00442 *
00443                   CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
00444      $                         'T', 'N', WORK( N+1 ), 1, ONE,
00445      $                         WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, K, K,
00446      $                         ZERO, ANORM, 'Q', A, LDA, IDUMMA, IINFO )
00447 *
00448                ELSE IF( ITYPE.EQ.9 ) THEN
00449 *
00450 *                 Positive definite, eigenvalues specified.
00451 *
00452                   CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
00453      $                         ANORM, K, K, 'Q', A, LDA, WORK( N+1 ),
00454      $                         IINFO )
00455 *
00456                ELSE IF( ITYPE.EQ.10 ) THEN
00457 *
00458 *                 Positive definite tridiagonal, eigenvalues specified.
00459 *
00460                   IF( N.GT.1 )
00461      $               K = MAX( 1, K )
00462                   CALL DLATMS( N, N, 'S', ISEED, 'P', WORK, IMODE, COND,
00463      $                         ANORM, 1, 1, 'Q', A( K, 1 ), LDA,
00464      $                         WORK( N+1 ), IINFO )
00465                   DO 90 I = 2, N
00466                      TEMP1 = ABS( A( K, I ) ) /
00467      $                       SQRT( ABS( A( K+1, I-1 )*A( K+1, I ) ) )
00468                      IF( TEMP1.GT.HALF ) THEN
00469                         A( K, I ) = HALF*SQRT( ABS( A( K+1,
00470      $                              I-1 )*A( K+1, I ) ) )
00471                      END IF
00472    90             CONTINUE
00473 *
00474                ELSE
00475 *
00476                   IINFO = 1
00477                END IF
00478 *
00479                IF( IINFO.NE.0 ) THEN
00480                   WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N,
00481      $               JTYPE, IOLDSD
00482                   INFO = ABS( IINFO )
00483                   RETURN
00484                END IF
00485 *
00486   100          CONTINUE
00487 *
00488 *              Call DSBTRD to compute S and U from upper triangle.
00489 *
00490                CALL DLACPY( ' ', K+1, N, A, LDA, WORK, LDA )
00491 *
00492                NTEST = 1
00493                CALL DSBTRD( 'V', 'U', N, K, WORK, LDA, SD, SE, U, LDU,
00494      $                      WORK( LDA*N+1 ), IINFO )
00495 *
00496                IF( IINFO.NE.0 ) THEN
00497                   WRITE( NOUNIT, FMT = 9999 )'DSBTRD(U)', IINFO, N,
00498      $               JTYPE, IOLDSD
00499                   INFO = ABS( IINFO )
00500                   IF( IINFO.LT.0 ) THEN
00501                      RETURN
00502                   ELSE
00503                      RESULT( 1 ) = ULPINV
00504                      GO TO 150
00505                   END IF
00506                END IF
00507 *
00508 *              Do tests 1 and 2
00509 *
00510                CALL DSBT21( 'Upper', N, K, 1, A, LDA, SD, SE, U, LDU,
00511      $                      WORK, RESULT( 1 ) )
00512 *
00513 *              Convert A from Upper-Triangle-Only storage to
00514 *              Lower-Triangle-Only storage.
00515 *
00516                DO 120 JC = 1, N
00517                   DO 110 JR = 0, MIN( K, N-JC )
00518                      A( JR+1, JC ) = A( K+1-JR, JC+JR )
00519   110             CONTINUE
00520   120          CONTINUE
00521                DO 140 JC = N + 1 - K, N
00522                   DO 130 JR = MIN( K, N-JC ) + 1, K
00523                      A( JR+1, JC ) = ZERO
00524   130             CONTINUE
00525   140          CONTINUE
00526 *
00527 *              Call DSBTRD to compute S and U from lower triangle
00528 *
00529                CALL DLACPY( ' ', K+1, N, A, LDA, WORK, LDA )
00530 *
00531                NTEST = 3
00532                CALL DSBTRD( 'V', 'L', N, K, WORK, LDA, SD, SE, U, LDU,
00533      $                      WORK( LDA*N+1 ), IINFO )
00534 *
00535                IF( IINFO.NE.0 ) THEN
00536                   WRITE( NOUNIT, FMT = 9999 )'DSBTRD(L)', IINFO, N,
00537      $               JTYPE, IOLDSD
00538                   INFO = ABS( IINFO )
00539                   IF( IINFO.LT.0 ) THEN
00540                      RETURN
00541                   ELSE
00542                      RESULT( 3 ) = ULPINV
00543                      GO TO 150
00544                   END IF
00545                END IF
00546                NTEST = 4
00547 *
00548 *              Do tests 3 and 4
00549 *
00550                CALL DSBT21( 'Lower', N, K, 1, A, LDA, SD, SE, U, LDU,
00551      $                      WORK, RESULT( 3 ) )
00552 *
00553 *              End of Loop -- Check for RESULT(j) > THRESH
00554 *
00555   150          CONTINUE
00556                NTESTT = NTESTT + NTEST
00557 *
00558 *              Print out tests which fail.
00559 *
00560                DO 160 JR = 1, NTEST
00561                   IF( RESULT( JR ).GE.THRESH ) THEN
00562 *
00563 *                    If this is the first test to fail,
00564 *                    print a header to the data file.
00565 *
00566                      IF( NERRS.EQ.0 ) THEN
00567                         WRITE( NOUNIT, FMT = 9998 )'DSB'
00568                         WRITE( NOUNIT, FMT = 9997 )
00569                         WRITE( NOUNIT, FMT = 9996 )
00570                         WRITE( NOUNIT, FMT = 9995 )'Symmetric'
00571                         WRITE( NOUNIT, FMT = 9994 )'orthogonal', '''',
00572      $                     'transpose', ( '''', J = 1, 4 )
00573                      END IF
00574                      NERRS = NERRS + 1
00575                      WRITE( NOUNIT, FMT = 9993 )N, K, IOLDSD, JTYPE,
00576      $                  JR, RESULT( JR )
00577                   END IF
00578   160          CONTINUE
00579 *
00580   170       CONTINUE
00581   180    CONTINUE
00582   190 CONTINUE
00583 *
00584 *     Summary
00585 *
00586       CALL DLASUM( 'DSB', NOUNIT, NERRS, NTESTT )
00587       RETURN
00588 *
00589  9999 FORMAT( ' DCHKSB: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00590      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
00591 *
00592  9998 FORMAT( / 1X, A3,
00593      $      ' -- Real Symmetric Banded Tridiagonal Reduction Routines' )
00594  9997 FORMAT( ' Matrix types (see DCHKSB for details): ' )
00595 *
00596  9996 FORMAT( / ' Special Matrices:',
00597      $      / '  1=Zero matrix.                        ',
00598      $      '  5=Diagonal: clustered entries.',
00599      $      / '  2=Identity matrix.                    ',
00600      $      '  6=Diagonal: large, evenly spaced.',
00601      $      / '  3=Diagonal: evenly spaced entries.    ',
00602      $      '  7=Diagonal: small, evenly spaced.',
00603      $      / '  4=Diagonal: geometr. spaced entries.' )
00604  9995 FORMAT( ' Dense ', A, ' Banded Matrices:',
00605      $      / '  8=Evenly spaced eigenvals.            ',
00606      $      ' 12=Small, evenly spaced eigenvals.',
00607      $      / '  9=Geometrically spaced eigenvals.     ',
00608      $      ' 13=Matrix with random O(1) entries.',
00609      $      / ' 10=Clustered eigenvalues.              ',
00610      $      ' 14=Matrix with large random entries.',
00611      $      / ' 11=Large, evenly spaced eigenvals.     ',
00612      $      ' 15=Matrix with small random entries.' )
00613 *
00614  9994 FORMAT( / ' Tests performed:   (S is Tridiag,  U is ', A, ',',
00615      $      / 20X, A, ' means ', A, '.', / ' UPLO=''U'':',
00616      $      / '  1= | A - U S U', A1, ' | / ( |A| n ulp )     ',
00617      $      '  2= | I - U U', A1, ' | / ( n ulp )', / ' UPLO=''L'':',
00618      $      / '  3= | A - U S U', A1, ' | / ( |A| n ulp )     ',
00619      $      '  4= | I - U U', A1, ' | / ( n ulp )' )
00620  9993 FORMAT( ' N=', I5, ', K=', I4, ', seed=', 4( I4, ',' ), ' type ',
00621      $      I2, ', test(', I2, ')=', G10.3 )
00622 *
00623 *     End of DCHKSB
00624 *
00625       END