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