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