 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGET38( 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 REAL RMAX( 3 ) 00013 * .. 00014 * 00015 * Purpose 00016 * ======= 00017 * 00018 * CGET38 tests CTRSEN, 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) REAL array, dimension (3) 00027 * Values of the largest test ratios. 00028 * RMAX(1) = largest residuals from CHST01 or comparing 00029 * different calls to CTRSEN 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 CGEHRD returns INFO nonzero on example i, LMAX(1)=i 00040 * If CHSEQR returns INFO nonzero on example i, LMAX(2)=i 00041 * If CTRSEN returns INFO nonzero on example i, LMAX(3)=i 00042 * 00043 * NINFO (output) INTEGER array, dimension (3) 00044 * NINFO(1) = No. of times CGEHRD returned INFO nonzero 00045 * NINFO(2) = No. of times CHSEQR returned INFO nonzero 00046 * NINFO(3) = No. of times CTRSEN 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 INTEGER LDT, LWORK 00058 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) ) 00059 REAL ZERO, ONE, TWO 00060 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 00061 REAL EPSIN 00062 PARAMETER ( EPSIN = 5.9605E-8 ) 00063 COMPLEX CZERO 00064 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00065 * .. 00066 * .. Local Scalars .. 00067 INTEGER I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM 00068 REAL BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN, 00069 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN, 00070 $ VMUL 00071 * .. 00072 * .. Local Arrays .. 00073 LOGICAL SELECT( LDT ) 00074 INTEGER IPNT( LDT ), ISELEC( LDT ) 00075 REAL RESULT( 2 ), RWORK( LDT ), VAL( 3 ), 00076 $ WSRT( LDT ) 00077 COMPLEX Q( LDT, LDT ), QSAV( LDT, LDT ), 00078 $ QTMP( LDT, LDT ), T( LDT, LDT ), 00079 $ TMP( LDT, LDT ), TSAV( LDT, LDT ), 00080 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ), 00081 $ WORK( LWORK ), WTMP( LDT ) 00082 * .. 00083 * .. External Functions .. 00084 REAL CLANGE, SLAMCH 00085 EXTERNAL CLANGE, SLAMCH 00086 * .. 00087 * .. External Subroutines .. 00088 EXTERNAL CGEHRD, CHSEQR, CHST01, CLACPY, CSSCAL, CTRSEN, 00089 $ CUNGHR, SLABAD 00090 * .. 00091 * .. Intrinsic Functions .. 00092 INTRINSIC AIMAG, MAX, REAL, SQRT 00093 * .. 00094 * .. Executable Statements .. 00095 * 00096 EPS = SLAMCH( 'P' ) 00097 SMLNUM = SLAMCH( 'S' ) / EPS 00098 BIGNUM = ONE / SMLNUM 00099 CALL SLABAD( SMLNUM, BIGNUM ) 00100 * 00101 * EPSIN = 2**(-24) = precision to which input data computed 00102 * 00103 EPS = MAX( EPS, EPSIN ) 00104 RMAX( 1 ) = ZERO 00105 RMAX( 2 ) = ZERO 00106 RMAX( 3 ) = ZERO 00107 LMAX( 1 ) = 0 00108 LMAX( 2 ) = 0 00109 LMAX( 3 ) = 0 00110 KNT = 0 00111 NINFO( 1 ) = 0 00112 NINFO( 2 ) = 0 00113 NINFO( 3 ) = 0 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, ISRT 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 = CLANGE( 'M', N, N, TMP, LDT, RWORK ) 00133 DO 200 ISCL = 1, 3 00134 * 00135 * Scale input matrix 00136 * 00137 KNT = KNT + 1 00138 CALL CLACPY( 'F', N, N, TMP, LDT, T, LDT ) 00139 VMUL = VAL( ISCL ) 00140 DO 30 I = 1, N 00141 CALL CSSCAL( N, VMUL, T( 1, I ), 1 ) 00142 30 CONTINUE 00143 IF( TNRM.EQ.ZERO ) 00144 $ VMUL = ONE 00145 CALL CLACPY( 'F', N, N, T, LDT, TSAV, LDT ) 00146 * 00147 * Compute Schur form 00148 * 00149 CALL CGEHRD( 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 200 00155 END IF 00156 * 00157 * Generate unitary matrix 00158 * 00159 CALL CLACPY( 'L', N, N, T, LDT, Q, LDT ) 00160 CALL CUNGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N, 00161 $ INFO ) 00162 * 00163 * Compute Schur form 00164 * 00165 DO 50 J = 1, N - 2 00166 DO 40 I = J + 2, N 00167 T( I, J ) = CZERO 00168 40 CONTINUE 00169 50 CONTINUE 00170 CALL CHSEQR( 'S', 'V', N, 1, N, T, LDT, W, Q, LDT, WORK, LWORK, 00171 $ INFO ) 00172 IF( INFO.NE.0 ) THEN 00173 LMAX( 2 ) = KNT 00174 NINFO( 2 ) = NINFO( 2 ) + 1 00175 GO TO 200 00176 END IF 00177 * 00178 * Sort, select eigenvalues 00179 * 00180 DO 60 I = 1, N 00181 IPNT( I ) = I 00182 SELECT( I ) = .FALSE. 00183 60 CONTINUE 00184 IF( ISRT.EQ.0 ) THEN 00185 DO 70 I = 1, N 00186 WSRT( I ) = REAL( W( I ) ) 00187 70 CONTINUE 00188 ELSE 00189 DO 80 I = 1, N 00190 WSRT( I ) = AIMAG( W( I ) ) 00191 80 CONTINUE 00192 END IF 00193 DO 100 I = 1, N - 1 00194 KMIN = I 00195 VMIN = WSRT( I ) 00196 DO 90 J = I + 1, N 00197 IF( WSRT( J ).LT.VMIN ) THEN 00198 KMIN = J 00199 VMIN = WSRT( J ) 00200 END IF 00201 90 CONTINUE 00202 WSRT( KMIN ) = WSRT( I ) 00203 WSRT( I ) = VMIN 00204 ITMP = IPNT( I ) 00205 IPNT( I ) = IPNT( KMIN ) 00206 IPNT( KMIN ) = ITMP 00207 100 CONTINUE 00208 DO 110 I = 1, NDIM 00209 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE. 00210 110 CONTINUE 00211 * 00212 * Compute condition numbers 00213 * 00214 CALL CLACPY( 'F', N, N, Q, LDT, QSAV, LDT ) 00215 CALL CLACPY( 'F', N, N, T, LDT, TSAV1, LDT ) 00216 CALL CTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WTMP, M, S, 00217 $ SEP, WORK, LWORK, INFO ) 00218 IF( INFO.NE.0 ) THEN 00219 LMAX( 3 ) = KNT 00220 NINFO( 3 ) = NINFO( 3 ) + 1 00221 GO TO 200 00222 END IF 00223 SEPTMP = SEP / VMUL 00224 STMP = S 00225 * 00226 * Compute residuals 00227 * 00228 CALL CHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK, 00229 $ RWORK, RESULT ) 00230 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) ) 00231 IF( VMAX.GT.RMAX( 1 ) ) THEN 00232 RMAX( 1 ) = VMAX 00233 IF( NINFO( 1 ).EQ.0 ) 00234 $ LMAX( 1 ) = KNT 00235 END IF 00236 * 00237 * Compare condition number for eigenvalue cluster 00238 * taking its condition number into account 00239 * 00240 V = MAX( TWO*REAL( N )*EPS*TNRM, SMLNUM ) 00241 IF( TNRM.EQ.ZERO ) 00242 $ V = ONE 00243 IF( V.GT.SEPTMP ) THEN 00244 TOL = ONE 00245 ELSE 00246 TOL = V / SEPTMP 00247 END IF 00248 IF( V.GT.SEPIN ) THEN 00249 TOLIN = ONE 00250 ELSE 00251 TOLIN = V / SEPIN 00252 END IF 00253 TOL = MAX( TOL, SMLNUM / EPS ) 00254 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00255 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN 00256 VMAX = ONE / EPS 00257 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN 00258 VMAX = ( SIN-TOLIN ) / ( STMP+TOL ) 00259 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN 00260 VMAX = ONE / EPS 00261 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN 00262 VMAX = ( STMP-TOL ) / ( SIN+TOLIN ) 00263 ELSE 00264 VMAX = ONE 00265 END IF 00266 IF( VMAX.GT.RMAX( 2 ) ) THEN 00267 RMAX( 2 ) = VMAX 00268 IF( NINFO( 2 ).EQ.0 ) 00269 $ LMAX( 2 ) = KNT 00270 END IF 00271 * 00272 * Compare condition numbers for invariant subspace 00273 * taking its condition number into account 00274 * 00275 IF( V.GT.SEPTMP*STMP ) THEN 00276 TOL = SEPTMP 00277 ELSE 00278 TOL = V / STMP 00279 END IF 00280 IF( V.GT.SEPIN*SIN ) THEN 00281 TOLIN = SEPIN 00282 ELSE 00283 TOLIN = V / SIN 00284 END IF 00285 TOL = MAX( TOL, SMLNUM / EPS ) 00286 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00287 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN 00288 VMAX = ONE / EPS 00289 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN 00290 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL ) 00291 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN 00292 VMAX = ONE / EPS 00293 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN 00294 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN ) 00295 ELSE 00296 VMAX = ONE 00297 END IF 00298 IF( VMAX.GT.RMAX( 2 ) ) THEN 00299 RMAX( 2 ) = VMAX 00300 IF( NINFO( 2 ).EQ.0 ) 00301 $ LMAX( 2 ) = KNT 00302 END IF 00303 * 00304 * Compare condition number for eigenvalue cluster 00305 * without taking its condition number into account 00306 * 00307 IF( SIN.LE.REAL( 2*N )*EPS .AND. STMP.LE.REAL( 2*N )*EPS ) THEN 00308 VMAX = ONE 00309 ELSE IF( EPS*SIN.GT.STMP ) THEN 00310 VMAX = ONE / EPS 00311 ELSE IF( SIN.GT.STMP ) THEN 00312 VMAX = SIN / STMP 00313 ELSE IF( SIN.LT.EPS*STMP ) THEN 00314 VMAX = ONE / EPS 00315 ELSE IF( SIN.LT.STMP ) THEN 00316 VMAX = STMP / SIN 00317 ELSE 00318 VMAX = ONE 00319 END IF 00320 IF( VMAX.GT.RMAX( 3 ) ) THEN 00321 RMAX( 3 ) = VMAX 00322 IF( NINFO( 3 ).EQ.0 ) 00323 $ LMAX( 3 ) = KNT 00324 END IF 00325 * 00326 * Compare condition numbers for invariant subspace 00327 * without taking its condition number into account 00328 * 00329 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN 00330 VMAX = ONE 00331 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN 00332 VMAX = ONE / EPS 00333 ELSE IF( SEPIN.GT.SEPTMP ) THEN 00334 VMAX = SEPIN / SEPTMP 00335 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN 00336 VMAX = ONE / EPS 00337 ELSE IF( SEPIN.LT.SEPTMP ) THEN 00338 VMAX = SEPTMP / SEPIN 00339 ELSE 00340 VMAX = ONE 00341 END IF 00342 IF( VMAX.GT.RMAX( 3 ) ) THEN 00343 RMAX( 3 ) = VMAX 00344 IF( NINFO( 3 ).EQ.0 ) 00345 $ LMAX( 3 ) = KNT 00346 END IF 00347 * 00348 * Compute eigenvalue condition number only and compare 00349 * Update Q 00350 * 00351 VMAX = ZERO 00352 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00353 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00354 SEPTMP = -ONE 00355 STMP = -ONE 00356 CALL CTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00357 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00358 IF( INFO.NE.0 ) THEN 00359 LMAX( 3 ) = KNT 00360 NINFO( 3 ) = NINFO( 3 ) + 1 00361 GO TO 200 00362 END IF 00363 IF( S.NE.STMP ) 00364 $ VMAX = ONE / EPS 00365 IF( -ONE.NE.SEPTMP ) 00366 $ VMAX = ONE / EPS 00367 DO 130 I = 1, N 00368 DO 120 J = 1, N 00369 IF( TTMP( I, J ).NE.T( I, J ) ) 00370 $ VMAX = ONE / EPS 00371 IF( QTMP( I, J ).NE.Q( I, J ) ) 00372 $ VMAX = ONE / EPS 00373 120 CONTINUE 00374 130 CONTINUE 00375 * 00376 * Compute invariant subspace condition number only and compare 00377 * Update Q 00378 * 00379 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00380 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00381 SEPTMP = -ONE 00382 STMP = -ONE 00383 CALL CTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00384 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00385 IF( INFO.NE.0 ) THEN 00386 LMAX( 3 ) = KNT 00387 NINFO( 3 ) = NINFO( 3 ) + 1 00388 GO TO 200 00389 END IF 00390 IF( -ONE.NE.STMP ) 00391 $ VMAX = ONE / EPS 00392 IF( SEP.NE.SEPTMP ) 00393 $ VMAX = ONE / EPS 00394 DO 150 I = 1, N 00395 DO 140 J = 1, N 00396 IF( TTMP( I, J ).NE.T( I, J ) ) 00397 $ VMAX = ONE / EPS 00398 IF( QTMP( I, J ).NE.Q( I, J ) ) 00399 $ VMAX = ONE / EPS 00400 140 CONTINUE 00401 150 CONTINUE 00402 * 00403 * Compute eigenvalue condition number only and compare 00404 * Do not update Q 00405 * 00406 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00407 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00408 SEPTMP = -ONE 00409 STMP = -ONE 00410 CALL CTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00411 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00412 IF( INFO.NE.0 ) THEN 00413 LMAX( 3 ) = KNT 00414 NINFO( 3 ) = NINFO( 3 ) + 1 00415 GO TO 200 00416 END IF 00417 IF( S.NE.STMP ) 00418 $ VMAX = ONE / EPS 00419 IF( -ONE.NE.SEPTMP ) 00420 $ VMAX = ONE / EPS 00421 DO 170 I = 1, N 00422 DO 160 J = 1, N 00423 IF( TTMP( I, J ).NE.T( I, J ) ) 00424 $ VMAX = ONE / EPS 00425 IF( QTMP( I, J ).NE.QSAV( I, J ) ) 00426 $ VMAX = ONE / EPS 00427 160 CONTINUE 00428 170 CONTINUE 00429 * 00430 * Compute invariant subspace condition number only and compare 00431 * Do not update Q 00432 * 00433 CALL CLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00434 CALL CLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00435 SEPTMP = -ONE 00436 STMP = -ONE 00437 CALL CTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00438 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00439 IF( INFO.NE.0 ) THEN 00440 LMAX( 3 ) = KNT 00441 NINFO( 3 ) = NINFO( 3 ) + 1 00442 GO TO 200 00443 END IF 00444 IF( -ONE.NE.STMP ) 00445 $ VMAX = ONE / EPS 00446 IF( SEP.NE.SEPTMP ) 00447 $ VMAX = ONE / EPS 00448 DO 190 I = 1, N 00449 DO 180 J = 1, N 00450 IF( TTMP( I, J ).NE.T( I, J ) ) 00451 $ VMAX = ONE / EPS 00452 IF( QTMP( I, J ).NE.QSAV( I, J ) ) 00453 $ VMAX = ONE / EPS 00454 180 CONTINUE 00455 190 CONTINUE 00456 IF( VMAX.GT.RMAX( 1 ) ) THEN 00457 RMAX( 1 ) = VMAX 00458 IF( NINFO( 1 ).EQ.0 ) 00459 $ LMAX( 1 ) = KNT 00460 END IF 00461 200 CONTINUE 00462 GO TO 10 00463 * 00464 * End of CGET38 00465 * 00466 END
1.7.3 
