![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGET38 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 DGET38( RMAX, LMAX, NINFO, KNT, NIN ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER KNT, NIN 00015 * .. 00016 * .. Array Arguments .. 00017 * INTEGER LMAX( 3 ), NINFO( 3 ) 00018 * DOUBLE PRECISION RMAX( 3 ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> DGET38 tests DTRSEN, a routine for estimating condition numbers of a 00028 *> cluster of eigenvalues and/or its associated right invariant subspace 00029 *> 00030 *> The test matrices are read from a file with logical unit number NIN. 00031 *> \endverbatim 00032 * 00033 * Arguments: 00034 * ========== 00035 * 00036 *> \param[out] RMAX 00037 *> \verbatim 00038 *> RMAX is DOUBLE PRECISION array, dimension (3) 00039 *> Values of the largest test ratios. 00040 *> RMAX(1) = largest residuals from DHST01 or comparing 00041 *> different calls to DTRSEN 00042 *> RMAX(2) = largest error in reciprocal condition 00043 *> numbers taking their conditioning into account 00044 *> RMAX(3) = largest error in reciprocal condition 00045 *> numbers not taking their conditioning into 00046 *> account (may be larger than RMAX(2)) 00047 *> \endverbatim 00048 *> 00049 *> \param[out] LMAX 00050 *> \verbatim 00051 *> LMAX is INTEGER array, dimension (3) 00052 *> LMAX(i) is example number where largest test ratio 00053 *> RMAX(i) is achieved. Also: 00054 *> If DGEHRD returns INFO nonzero on example i, LMAX(1)=i 00055 *> If DHSEQR returns INFO nonzero on example i, LMAX(2)=i 00056 *> If DTRSEN returns INFO nonzero on example i, LMAX(3)=i 00057 *> \endverbatim 00058 *> 00059 *> \param[out] NINFO 00060 *> \verbatim 00061 *> NINFO is INTEGER array, dimension (3) 00062 *> NINFO(1) = No. of times DGEHRD returned INFO nonzero 00063 *> NINFO(2) = No. of times DHSEQR returned INFO nonzero 00064 *> NINFO(3) = No. of times DTRSEN returned INFO nonzero 00065 *> \endverbatim 00066 *> 00067 *> \param[out] KNT 00068 *> \verbatim 00069 *> KNT is INTEGER 00070 *> Total number of examples tested. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] NIN 00074 *> \verbatim 00075 *> NIN is INTEGER 00076 *> Input logical unit number. 00077 *> \endverbatim 00078 * 00079 * Authors: 00080 * ======== 00081 * 00082 *> \author Univ. of Tennessee 00083 *> \author Univ. of California Berkeley 00084 *> \author Univ. of Colorado Denver 00085 *> \author NAG Ltd. 00086 * 00087 *> \date November 2011 00088 * 00089 *> \ingroup double_eig 00090 * 00091 * ===================================================================== 00092 SUBROUTINE DGET38( RMAX, LMAX, NINFO, KNT, NIN ) 00093 * 00094 * -- LAPACK test routine (version 3.4.0) -- 00095 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00096 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00097 * November 2011 00098 * 00099 * .. Scalar Arguments .. 00100 INTEGER KNT, NIN 00101 * .. 00102 * .. Array Arguments .. 00103 INTEGER LMAX( 3 ), NINFO( 3 ) 00104 DOUBLE PRECISION RMAX( 3 ) 00105 * .. 00106 * 00107 * ===================================================================== 00108 * 00109 * .. Parameters .. 00110 DOUBLE PRECISION ZERO, ONE, TWO 00111 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00112 DOUBLE PRECISION EPSIN 00113 PARAMETER ( EPSIN = 5.9605D-8 ) 00114 INTEGER LDT, LWORK 00115 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) ) 00116 INTEGER LIWORK 00117 PARAMETER ( LIWORK = LDT*LDT ) 00118 * .. 00119 * .. Local Scalars .. 00120 INTEGER I, INFO, ISCL, ITMP, J, KMIN, M, N, NDIM 00121 DOUBLE PRECISION BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN, 00122 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VIMIN, VMAX, 00123 $ VMUL, VRMIN 00124 * .. 00125 * .. Local Arrays .. 00126 LOGICAL SELECT( LDT ) 00127 INTEGER IPNT( LDT ), ISELEC( LDT ), IWORK( LIWORK ) 00128 DOUBLE PRECISION Q( LDT, LDT ), QSAV( LDT, LDT ), 00129 $ QTMP( LDT, LDT ), RESULT( 2 ), T( LDT, LDT ), 00130 $ TMP( LDT, LDT ), TSAV( LDT, LDT ), 00131 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), VAL( 3 ), 00132 $ WI( LDT ), WITMP( LDT ), WORK( LWORK ), 00133 $ WR( LDT ), WRTMP( LDT ) 00134 * .. 00135 * .. External Functions .. 00136 DOUBLE PRECISION DLAMCH, DLANGE 00137 EXTERNAL DLAMCH, DLANGE 00138 * .. 00139 * .. External Subroutines .. 00140 EXTERNAL DCOPY, DGEHRD, DHSEQR, DHST01, DLABAD, DLACPY, 00141 $ DORGHR, DSCAL, DTRSEN 00142 * .. 00143 * .. Intrinsic Functions .. 00144 INTRINSIC DBLE, MAX, SQRT 00145 * .. 00146 * .. Executable Statements .. 00147 * 00148 EPS = DLAMCH( 'P' ) 00149 SMLNUM = DLAMCH( 'S' ) / EPS 00150 BIGNUM = ONE / SMLNUM 00151 CALL DLABAD( SMLNUM, BIGNUM ) 00152 * 00153 * EPSIN = 2**(-24) = precision to which input data computed 00154 * 00155 EPS = MAX( EPS, EPSIN ) 00156 RMAX( 1 ) = ZERO 00157 RMAX( 2 ) = ZERO 00158 RMAX( 3 ) = ZERO 00159 LMAX( 1 ) = 0 00160 LMAX( 2 ) = 0 00161 LMAX( 3 ) = 0 00162 KNT = 0 00163 NINFO( 1 ) = 0 00164 NINFO( 2 ) = 0 00165 NINFO( 3 ) = 0 00166 * 00167 VAL( 1 ) = SQRT( SMLNUM ) 00168 VAL( 2 ) = ONE 00169 VAL( 3 ) = SQRT( SQRT( BIGNUM ) ) 00170 * 00171 * Read input data until N=0. Assume input eigenvalues are sorted 00172 * lexicographically (increasing by real part, then decreasing by 00173 * imaginary part) 00174 * 00175 10 CONTINUE 00176 READ( NIN, FMT = * )N, NDIM 00177 IF( N.EQ.0 ) 00178 $ RETURN 00179 READ( NIN, FMT = * )( ISELEC( I ), I = 1, NDIM ) 00180 DO 20 I = 1, N 00181 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N ) 00182 20 CONTINUE 00183 READ( NIN, FMT = * )SIN, SEPIN 00184 * 00185 TNRM = DLANGE( 'M', N, N, TMP, LDT, WORK ) 00186 DO 160 ISCL = 1, 3 00187 * 00188 * Scale input matrix 00189 * 00190 KNT = KNT + 1 00191 CALL DLACPY( 'F', N, N, TMP, LDT, T, LDT ) 00192 VMUL = VAL( ISCL ) 00193 DO 30 I = 1, N 00194 CALL DSCAL( N, VMUL, T( 1, I ), 1 ) 00195 30 CONTINUE 00196 IF( TNRM.EQ.ZERO ) 00197 $ VMUL = ONE 00198 CALL DLACPY( 'F', N, N, T, LDT, TSAV, LDT ) 00199 * 00200 * Compute Schur form 00201 * 00202 CALL DGEHRD( N, 1, N, T, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N, 00203 $ INFO ) 00204 IF( INFO.NE.0 ) THEN 00205 LMAX( 1 ) = KNT 00206 NINFO( 1 ) = NINFO( 1 ) + 1 00207 GO TO 160 00208 END IF 00209 * 00210 * Generate orthogonal matrix 00211 * 00212 CALL DLACPY( 'L', N, N, T, LDT, Q, LDT ) 00213 CALL DORGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N, 00214 $ INFO ) 00215 * 00216 * Compute Schur form 00217 * 00218 CALL DHSEQR( 'S', 'V', N, 1, N, T, LDT, WR, WI, Q, LDT, WORK, 00219 $ LWORK, INFO ) 00220 IF( INFO.NE.0 ) THEN 00221 LMAX( 2 ) = KNT 00222 NINFO( 2 ) = NINFO( 2 ) + 1 00223 GO TO 160 00224 END IF 00225 * 00226 * Sort, select eigenvalues 00227 * 00228 DO 40 I = 1, N 00229 IPNT( I ) = I 00230 SELECT( I ) = .FALSE. 00231 40 CONTINUE 00232 CALL DCOPY( N, WR, 1, WRTMP, 1 ) 00233 CALL DCOPY( N, WI, 1, WITMP, 1 ) 00234 DO 60 I = 1, N - 1 00235 KMIN = I 00236 VRMIN = WRTMP( I ) 00237 VIMIN = WITMP( I ) 00238 DO 50 J = I + 1, N 00239 IF( WRTMP( J ).LT.VRMIN ) THEN 00240 KMIN = J 00241 VRMIN = WRTMP( J ) 00242 VIMIN = WITMP( J ) 00243 END IF 00244 50 CONTINUE 00245 WRTMP( KMIN ) = WRTMP( I ) 00246 WITMP( KMIN ) = WITMP( I ) 00247 WRTMP( I ) = VRMIN 00248 WITMP( I ) = VIMIN 00249 ITMP = IPNT( I ) 00250 IPNT( I ) = IPNT( KMIN ) 00251 IPNT( KMIN ) = ITMP 00252 60 CONTINUE 00253 DO 70 I = 1, NDIM 00254 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE. 00255 70 CONTINUE 00256 * 00257 * Compute condition numbers 00258 * 00259 CALL DLACPY( 'F', N, N, Q, LDT, QSAV, LDT ) 00260 CALL DLACPY( 'F', N, N, T, LDT, TSAV1, LDT ) 00261 CALL DTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WRTMP, WITMP, 00262 $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO ) 00263 IF( INFO.NE.0 ) THEN 00264 LMAX( 3 ) = KNT 00265 NINFO( 3 ) = NINFO( 3 ) + 1 00266 GO TO 160 00267 END IF 00268 SEPTMP = SEP / VMUL 00269 STMP = S 00270 * 00271 * Compute residuals 00272 * 00273 CALL DHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK, 00274 $ RESULT ) 00275 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) ) 00276 IF( VMAX.GT.RMAX( 1 ) ) THEN 00277 RMAX( 1 ) = VMAX 00278 IF( NINFO( 1 ).EQ.0 ) 00279 $ LMAX( 1 ) = KNT 00280 END IF 00281 * 00282 * Compare condition number for eigenvalue cluster 00283 * taking its condition number into account 00284 * 00285 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM ) 00286 IF( TNRM.EQ.ZERO ) 00287 $ V = ONE 00288 IF( V.GT.SEPTMP ) THEN 00289 TOL = ONE 00290 ELSE 00291 TOL = V / SEPTMP 00292 END IF 00293 IF( V.GT.SEPIN ) THEN 00294 TOLIN = ONE 00295 ELSE 00296 TOLIN = V / SEPIN 00297 END IF 00298 TOL = MAX( TOL, SMLNUM / EPS ) 00299 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00300 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN 00301 VMAX = ONE / EPS 00302 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN 00303 VMAX = ( SIN-TOLIN ) / ( STMP+TOL ) 00304 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN 00305 VMAX = ONE / EPS 00306 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN 00307 VMAX = ( STMP-TOL ) / ( SIN+TOLIN ) 00308 ELSE 00309 VMAX = ONE 00310 END IF 00311 IF( VMAX.GT.RMAX( 2 ) ) THEN 00312 RMAX( 2 ) = VMAX 00313 IF( NINFO( 2 ).EQ.0 ) 00314 $ LMAX( 2 ) = KNT 00315 END IF 00316 * 00317 * Compare condition numbers for invariant subspace 00318 * taking its condition number into account 00319 * 00320 IF( V.GT.SEPTMP*STMP ) THEN 00321 TOL = SEPTMP 00322 ELSE 00323 TOL = V / STMP 00324 END IF 00325 IF( V.GT.SEPIN*SIN ) THEN 00326 TOLIN = SEPIN 00327 ELSE 00328 TOLIN = V / SIN 00329 END IF 00330 TOL = MAX( TOL, SMLNUM / EPS ) 00331 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00332 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN 00333 VMAX = ONE / EPS 00334 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN 00335 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL ) 00336 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN 00337 VMAX = ONE / EPS 00338 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN 00339 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN ) 00340 ELSE 00341 VMAX = ONE 00342 END IF 00343 IF( VMAX.GT.RMAX( 2 ) ) THEN 00344 RMAX( 2 ) = VMAX 00345 IF( NINFO( 2 ).EQ.0 ) 00346 $ LMAX( 2 ) = KNT 00347 END IF 00348 * 00349 * Compare condition number for eigenvalue cluster 00350 * without taking its condition number into account 00351 * 00352 IF( SIN.LE.DBLE( 2*N )*EPS .AND. STMP.LE.DBLE( 2*N )*EPS ) THEN 00353 VMAX = ONE 00354 ELSE IF( EPS*SIN.GT.STMP ) THEN 00355 VMAX = ONE / EPS 00356 ELSE IF( SIN.GT.STMP ) THEN 00357 VMAX = SIN / STMP 00358 ELSE IF( SIN.LT.EPS*STMP ) THEN 00359 VMAX = ONE / EPS 00360 ELSE IF( SIN.LT.STMP ) THEN 00361 VMAX = STMP / SIN 00362 ELSE 00363 VMAX = ONE 00364 END IF 00365 IF( VMAX.GT.RMAX( 3 ) ) THEN 00366 RMAX( 3 ) = VMAX 00367 IF( NINFO( 3 ).EQ.0 ) 00368 $ LMAX( 3 ) = KNT 00369 END IF 00370 * 00371 * Compare condition numbers for invariant subspace 00372 * without taking its condition number into account 00373 * 00374 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN 00375 VMAX = ONE 00376 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN 00377 VMAX = ONE / EPS 00378 ELSE IF( SEPIN.GT.SEPTMP ) THEN 00379 VMAX = SEPIN / SEPTMP 00380 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN 00381 VMAX = ONE / EPS 00382 ELSE IF( SEPIN.LT.SEPTMP ) THEN 00383 VMAX = SEPTMP / SEPIN 00384 ELSE 00385 VMAX = ONE 00386 END IF 00387 IF( VMAX.GT.RMAX( 3 ) ) THEN 00388 RMAX( 3 ) = VMAX 00389 IF( NINFO( 3 ).EQ.0 ) 00390 $ LMAX( 3 ) = KNT 00391 END IF 00392 * 00393 * Compute eigenvalue condition number only and compare 00394 * Update Q 00395 * 00396 VMAX = ZERO 00397 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00398 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00399 SEPTMP = -ONE 00400 STMP = -ONE 00401 CALL DTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP, 00402 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK, 00403 $ LIWORK, INFO ) 00404 IF( INFO.NE.0 ) THEN 00405 LMAX( 3 ) = KNT 00406 NINFO( 3 ) = NINFO( 3 ) + 1 00407 GO TO 160 00408 END IF 00409 IF( S.NE.STMP ) 00410 $ VMAX = ONE / EPS 00411 IF( -ONE.NE.SEPTMP ) 00412 $ VMAX = ONE / EPS 00413 DO 90 I = 1, N 00414 DO 80 J = 1, N 00415 IF( TTMP( I, J ).NE.T( I, J ) ) 00416 $ VMAX = ONE / EPS 00417 IF( QTMP( I, J ).NE.Q( I, J ) ) 00418 $ VMAX = ONE / EPS 00419 80 CONTINUE 00420 90 CONTINUE 00421 * 00422 * Compute invariant subspace condition number only and compare 00423 * Update Q 00424 * 00425 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00426 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00427 SEPTMP = -ONE 00428 STMP = -ONE 00429 CALL DTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP, 00430 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK, 00431 $ LIWORK, INFO ) 00432 IF( INFO.NE.0 ) THEN 00433 LMAX( 3 ) = KNT 00434 NINFO( 3 ) = NINFO( 3 ) + 1 00435 GO TO 160 00436 END IF 00437 IF( -ONE.NE.STMP ) 00438 $ VMAX = ONE / EPS 00439 IF( SEP.NE.SEPTMP ) 00440 $ VMAX = ONE / EPS 00441 DO 110 I = 1, N 00442 DO 100 J = 1, N 00443 IF( TTMP( I, J ).NE.T( I, J ) ) 00444 $ VMAX = ONE / EPS 00445 IF( QTMP( I, J ).NE.Q( I, J ) ) 00446 $ VMAX = ONE / EPS 00447 100 CONTINUE 00448 110 CONTINUE 00449 * 00450 * Compute eigenvalue condition number only and compare 00451 * Do not update Q 00452 * 00453 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00454 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00455 SEPTMP = -ONE 00456 STMP = -ONE 00457 CALL DTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP, 00458 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK, 00459 $ LIWORK, INFO ) 00460 IF( INFO.NE.0 ) THEN 00461 LMAX( 3 ) = KNT 00462 NINFO( 3 ) = NINFO( 3 ) + 1 00463 GO TO 160 00464 END IF 00465 IF( S.NE.STMP ) 00466 $ VMAX = ONE / EPS 00467 IF( -ONE.NE.SEPTMP ) 00468 $ VMAX = ONE / EPS 00469 DO 130 I = 1, N 00470 DO 120 J = 1, N 00471 IF( TTMP( I, J ).NE.T( I, J ) ) 00472 $ VMAX = ONE / EPS 00473 IF( QTMP( I, J ).NE.QSAV( I, J ) ) 00474 $ VMAX = ONE / EPS 00475 120 CONTINUE 00476 130 CONTINUE 00477 * 00478 * Compute invariant subspace condition number only and compare 00479 * Do not update Q 00480 * 00481 CALL DLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00482 CALL DLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00483 SEPTMP = -ONE 00484 STMP = -ONE 00485 CALL DTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP, 00486 $ WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK, 00487 $ LIWORK, INFO ) 00488 IF( INFO.NE.0 ) THEN 00489 LMAX( 3 ) = KNT 00490 NINFO( 3 ) = NINFO( 3 ) + 1 00491 GO TO 160 00492 END IF 00493 IF( -ONE.NE.STMP ) 00494 $ VMAX = ONE / EPS 00495 IF( SEP.NE.SEPTMP ) 00496 $ VMAX = ONE / EPS 00497 DO 150 I = 1, N 00498 DO 140 J = 1, N 00499 IF( TTMP( I, J ).NE.T( I, J ) ) 00500 $ VMAX = ONE / EPS 00501 IF( QTMP( I, J ).NE.QSAV( I, J ) ) 00502 $ VMAX = ONE / EPS 00503 140 CONTINUE 00504 150 CONTINUE 00505 IF( VMAX.GT.RMAX( 1 ) ) THEN 00506 RMAX( 1 ) = VMAX 00507 IF( NINFO( 1 ).EQ.0 ) 00508 $ LMAX( 1 ) = KNT 00509 END IF 00510 160 CONTINUE 00511 GO TO 10 00512 * 00513 * End of DGET38 00514 * 00515 END