![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DCHKSB 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 DCHKSB( 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 * DOUBLE PRECISION THRESH 00019 * .. 00020 * .. Array Arguments .. 00021 * LOGICAL DOTYPE( * ) 00022 * INTEGER ISEED( 4 ), KK( * ), NN( * ) 00023 * DOUBLE PRECISION A( LDA, * ), RESULT( * ), SD( * ), SE( * ), 00024 * $ U( LDU, * ), WORK( * ) 00025 * .. 00026 * 00027 * 00028 *> \par Purpose: 00029 * ============= 00030 *> 00031 *> \verbatim 00032 *> 00033 *> DCHKSB tests the reduction of a symmetric band matrix to tridiagonal 00034 *> form, used with the symmetric eigenvalue problem. 00035 *> 00036 *> DSBTRD factors a symmetric band matrix A as U S U' , where ' means 00037 *> transpose, S is symmetric tridiagonal, and U is orthogonal. 00038 *> DSBTRD can use either just the lower or just the upper triangle 00039 *> of A; DCHKSB checks both cases. 00040 *> 00041 *> When DCHKSB 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 DSBTRD with 00049 *> UPLO='U' 00050 *> 00051 *> (2) | I - UU' | / ( n ulp ) 00052 *> 00053 *> (3) | A - V S V' | / ( |A| n ulp ) computed by DSBTRD 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 *> DCHKSB 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 *> DCHKSB 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, DCHKSB 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 DCHKSB to continue the same random number 00164 *> sequence. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] THRESH 00168 *> \verbatim 00169 *> THRESH is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (max(NN)) 00203 *> Used to hold the diagonal of the tridiagonal matrix computed 00204 *> by DSBTRD. 00205 *> \endverbatim 00206 *> 00207 *> \param[out] SE 00208 *> \verbatim 00209 *> SE is DOUBLE PRECISION array, dimension (max(NN)) 00210 *> Used to hold the off-diagonal of the tridiagonal matrix 00211 *> computed by DSBTRD. 00212 *> \endverbatim 00213 *> 00214 *> \param[out] U 00215 *> \verbatim 00216 *> U is DOUBLE PRECISION array, dimension (LDU, max(NN)) 00217 *> Used to hold the orthogonal matrix computed by DSBTRD. 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 double_eig 00290 * 00291 * ===================================================================== 00292 SUBROUTINE DCHKSB( 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 DOUBLE PRECISION THRESH 00305 * .. 00306 * .. Array Arguments .. 00307 LOGICAL DOTYPE( * ) 00308 INTEGER ISEED( 4 ), KK( * ), NN( * ) 00309 DOUBLE PRECISION A( LDA, * ), RESULT( * ), SD( * ), SE( * ), 00310 $ U( LDU, * ), WORK( * ) 00311 * .. 00312 * 00313 * ===================================================================== 00314 * 00315 * .. Parameters .. 00316 DOUBLE PRECISION ZERO, ONE, TWO, TEN 00317 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00318 $ TEN = 10.0D0 ) 00319 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00338 EXTERNAL DLAMCH 00339 * .. 00340 * .. External Subroutines .. 00341 EXTERNAL DLACPY, DLASET, DLASUM, DLATMR, DLATMS, DSBT21, 00342 $ DSBTRD, XERBLA 00343 * .. 00344 * .. Intrinsic Functions .. 00345 INTRINSIC ABS, DBLE, MAX, MIN, 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( 'DCHKSB', -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 = DLAMCH( 'Safe minimum' ) 00413 OVFL = ONE / UNFL 00414 ULP = DLAMCH( 'Epsilon' )*DLAMCH( '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 / DBLE( 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 DLASET( '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 DLATMS( 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 DLATMS( 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 DLATMR( 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 DLATMR( 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 DLATMS( 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 DLATMS( 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 DSBTRD to compute S and U from upper triangle. 00590 * 00591 CALL DLACPY( ' ', K+1, N, A, LDA, WORK, LDA ) 00592 * 00593 NTEST = 1 00594 CALL DSBTRD( '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 )'DSBTRD(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 DSBT21( '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 DSBTRD to compute S and U from lower triangle 00629 * 00630 CALL DLACPY( ' ', K+1, N, A, LDA, WORK, LDA ) 00631 * 00632 NTEST = 3 00633 CALL DSBTRD( '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 )'DSBTRD(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 DSBT21( '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 )'DSB' 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 DLASUM( 'DSB', NOUNIT, NERRS, NTESTT ) 00688 RETURN 00689 * 00690 9999 FORMAT( ' DCHKSB: ', 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 DCHKSB 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 DCHKSB 00725 * 00726 END