![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGET38 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 ZGET38( 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 *> ZGET38 tests ZTRSEN, 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 ZHST01 or comparing 00041 *> different calls to ZTRSEN 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 ZGEHRD returns INFO nonzero on example i, LMAX(1)=i 00055 *> If ZHSEQR returns INFO nonzero on example i, LMAX(2)=i 00056 *> If ZTRSEN 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 ZGEHRD returned INFO nonzero 00063 *> NINFO(2) = No. of times ZHSEQR returned INFO nonzero 00064 *> NINFO(3) = No. of times ZTRSEN 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 complex16_eig 00090 * 00091 * ===================================================================== 00092 SUBROUTINE ZGET38( 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 INTEGER LDT, LWORK 00111 PARAMETER ( LDT = 20, LWORK = 2*LDT*( 10+LDT ) ) 00112 DOUBLE PRECISION ZERO, ONE, TWO 00113 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00114 DOUBLE PRECISION EPSIN 00115 PARAMETER ( EPSIN = 5.9605D-8 ) 00116 COMPLEX*16 CZERO 00117 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 00118 * .. 00119 * .. Local Scalars .. 00120 INTEGER I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM 00121 DOUBLE PRECISION BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN, 00122 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN, 00123 $ VMUL 00124 * .. 00125 * .. Local Arrays .. 00126 LOGICAL SELECT( LDT ) 00127 INTEGER IPNT( LDT ), ISELEC( LDT ) 00128 DOUBLE PRECISION RESULT( 2 ), RWORK( LDT ), VAL( 3 ), 00129 $ WSRT( LDT ) 00130 COMPLEX*16 Q( LDT, LDT ), QSAV( LDT, LDT ), 00131 $ QTMP( LDT, LDT ), T( LDT, LDT ), 00132 $ TMP( LDT, LDT ), TSAV( LDT, LDT ), 00133 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ), 00134 $ WORK( LWORK ), WTMP( LDT ) 00135 * .. 00136 * .. External Functions .. 00137 DOUBLE PRECISION DLAMCH, ZLANGE 00138 EXTERNAL DLAMCH, ZLANGE 00139 * .. 00140 * .. External Subroutines .. 00141 EXTERNAL DLABAD, ZDSCAL, ZGEHRD, ZHSEQR, ZHST01, ZLACPY, 00142 $ ZTRSEN, ZUNGHR 00143 * .. 00144 * .. Intrinsic Functions .. 00145 INTRINSIC DBLE, DIMAG, MAX, SQRT 00146 * .. 00147 * .. Executable Statements .. 00148 * 00149 EPS = DLAMCH( 'P' ) 00150 SMLNUM = DLAMCH( 'S' ) / EPS 00151 BIGNUM = ONE / SMLNUM 00152 CALL DLABAD( SMLNUM, BIGNUM ) 00153 * 00154 * EPSIN = 2**(-24) = precision to which input data computed 00155 * 00156 EPS = MAX( EPS, EPSIN ) 00157 RMAX( 1 ) = ZERO 00158 RMAX( 2 ) = ZERO 00159 RMAX( 3 ) = ZERO 00160 LMAX( 1 ) = 0 00161 LMAX( 2 ) = 0 00162 LMAX( 3 ) = 0 00163 KNT = 0 00164 NINFO( 1 ) = 0 00165 NINFO( 2 ) = 0 00166 NINFO( 3 ) = 0 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, ISRT 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 = ZLANGE( 'M', N, N, TMP, LDT, RWORK ) 00186 DO 200 ISCL = 1, 3 00187 * 00188 * Scale input matrix 00189 * 00190 KNT = KNT + 1 00191 CALL ZLACPY( 'F', N, N, TMP, LDT, T, LDT ) 00192 VMUL = VAL( ISCL ) 00193 DO 30 I = 1, N 00194 CALL ZDSCAL( N, VMUL, T( 1, I ), 1 ) 00195 30 CONTINUE 00196 IF( TNRM.EQ.ZERO ) 00197 $ VMUL = ONE 00198 CALL ZLACPY( 'F', N, N, T, LDT, TSAV, LDT ) 00199 * 00200 * Compute Schur form 00201 * 00202 CALL ZGEHRD( 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 200 00208 END IF 00209 * 00210 * Generate unitary matrix 00211 * 00212 CALL ZLACPY( 'L', N, N, T, LDT, Q, LDT ) 00213 CALL ZUNGHR( N, 1, N, Q, LDT, WORK( 1 ), WORK( N+1 ), LWORK-N, 00214 $ INFO ) 00215 * 00216 * Compute Schur form 00217 * 00218 DO 50 J = 1, N - 2 00219 DO 40 I = J + 2, N 00220 T( I, J ) = CZERO 00221 40 CONTINUE 00222 50 CONTINUE 00223 CALL ZHSEQR( 'S', 'V', N, 1, N, T, LDT, W, Q, LDT, WORK, LWORK, 00224 $ INFO ) 00225 IF( INFO.NE.0 ) THEN 00226 LMAX( 2 ) = KNT 00227 NINFO( 2 ) = NINFO( 2 ) + 1 00228 GO TO 200 00229 END IF 00230 * 00231 * Sort, select eigenvalues 00232 * 00233 DO 60 I = 1, N 00234 IPNT( I ) = I 00235 SELECT( I ) = .FALSE. 00236 60 CONTINUE 00237 IF( ISRT.EQ.0 ) THEN 00238 DO 70 I = 1, N 00239 WSRT( I ) = DBLE( W( I ) ) 00240 70 CONTINUE 00241 ELSE 00242 DO 80 I = 1, N 00243 WSRT( I ) = DIMAG( W( I ) ) 00244 80 CONTINUE 00245 END IF 00246 DO 100 I = 1, N - 1 00247 KMIN = I 00248 VMIN = WSRT( I ) 00249 DO 90 J = I + 1, N 00250 IF( WSRT( J ).LT.VMIN ) THEN 00251 KMIN = J 00252 VMIN = WSRT( J ) 00253 END IF 00254 90 CONTINUE 00255 WSRT( KMIN ) = WSRT( I ) 00256 WSRT( I ) = VMIN 00257 ITMP = IPNT( I ) 00258 IPNT( I ) = IPNT( KMIN ) 00259 IPNT( KMIN ) = ITMP 00260 100 CONTINUE 00261 DO 110 I = 1, NDIM 00262 SELECT( IPNT( ISELEC( I ) ) ) = .TRUE. 00263 110 CONTINUE 00264 * 00265 * Compute condition numbers 00266 * 00267 CALL ZLACPY( 'F', N, N, Q, LDT, QSAV, LDT ) 00268 CALL ZLACPY( 'F', N, N, T, LDT, TSAV1, LDT ) 00269 CALL ZTRSEN( 'B', 'V', SELECT, N, T, LDT, Q, LDT, WTMP, M, S, 00270 $ SEP, WORK, LWORK, INFO ) 00271 IF( INFO.NE.0 ) THEN 00272 LMAX( 3 ) = KNT 00273 NINFO( 3 ) = NINFO( 3 ) + 1 00274 GO TO 200 00275 END IF 00276 SEPTMP = SEP / VMUL 00277 STMP = S 00278 * 00279 * Compute residuals 00280 * 00281 CALL ZHST01( N, 1, N, TSAV, LDT, T, LDT, Q, LDT, WORK, LWORK, 00282 $ RWORK, RESULT ) 00283 VMAX = MAX( RESULT( 1 ), RESULT( 2 ) ) 00284 IF( VMAX.GT.RMAX( 1 ) ) THEN 00285 RMAX( 1 ) = VMAX 00286 IF( NINFO( 1 ).EQ.0 ) 00287 $ LMAX( 1 ) = KNT 00288 END IF 00289 * 00290 * Compare condition number for eigenvalue cluster 00291 * taking its condition number into account 00292 * 00293 V = MAX( TWO*DBLE( N )*EPS*TNRM, SMLNUM ) 00294 IF( TNRM.EQ.ZERO ) 00295 $ V = ONE 00296 IF( V.GT.SEPTMP ) THEN 00297 TOL = ONE 00298 ELSE 00299 TOL = V / SEPTMP 00300 END IF 00301 IF( V.GT.SEPIN ) THEN 00302 TOLIN = ONE 00303 ELSE 00304 TOLIN = V / SEPIN 00305 END IF 00306 TOL = MAX( TOL, SMLNUM / EPS ) 00307 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00308 IF( EPS*( SIN-TOLIN ).GT.STMP+TOL ) THEN 00309 VMAX = ONE / EPS 00310 ELSE IF( SIN-TOLIN.GT.STMP+TOL ) THEN 00311 VMAX = ( SIN-TOLIN ) / ( STMP+TOL ) 00312 ELSE IF( SIN+TOLIN.LT.EPS*( STMP-TOL ) ) THEN 00313 VMAX = ONE / EPS 00314 ELSE IF( SIN+TOLIN.LT.STMP-TOL ) THEN 00315 VMAX = ( STMP-TOL ) / ( SIN+TOLIN ) 00316 ELSE 00317 VMAX = ONE 00318 END IF 00319 IF( VMAX.GT.RMAX( 2 ) ) THEN 00320 RMAX( 2 ) = VMAX 00321 IF( NINFO( 2 ).EQ.0 ) 00322 $ LMAX( 2 ) = KNT 00323 END IF 00324 * 00325 * Compare condition numbers for invariant subspace 00326 * taking its condition number into account 00327 * 00328 IF( V.GT.SEPTMP*STMP ) THEN 00329 TOL = SEPTMP 00330 ELSE 00331 TOL = V / STMP 00332 END IF 00333 IF( V.GT.SEPIN*SIN ) THEN 00334 TOLIN = SEPIN 00335 ELSE 00336 TOLIN = V / SIN 00337 END IF 00338 TOL = MAX( TOL, SMLNUM / EPS ) 00339 TOLIN = MAX( TOLIN, SMLNUM / EPS ) 00340 IF( EPS*( SEPIN-TOLIN ).GT.SEPTMP+TOL ) THEN 00341 VMAX = ONE / EPS 00342 ELSE IF( SEPIN-TOLIN.GT.SEPTMP+TOL ) THEN 00343 VMAX = ( SEPIN-TOLIN ) / ( SEPTMP+TOL ) 00344 ELSE IF( SEPIN+TOLIN.LT.EPS*( SEPTMP-TOL ) ) THEN 00345 VMAX = ONE / EPS 00346 ELSE IF( SEPIN+TOLIN.LT.SEPTMP-TOL ) THEN 00347 VMAX = ( SEPTMP-TOL ) / ( SEPIN+TOLIN ) 00348 ELSE 00349 VMAX = ONE 00350 END IF 00351 IF( VMAX.GT.RMAX( 2 ) ) THEN 00352 RMAX( 2 ) = VMAX 00353 IF( NINFO( 2 ).EQ.0 ) 00354 $ LMAX( 2 ) = KNT 00355 END IF 00356 * 00357 * Compare condition number for eigenvalue cluster 00358 * without taking its condition number into account 00359 * 00360 IF( SIN.LE.DBLE( 2*N )*EPS .AND. STMP.LE.DBLE( 2*N )*EPS ) THEN 00361 VMAX = ONE 00362 ELSE IF( EPS*SIN.GT.STMP ) THEN 00363 VMAX = ONE / EPS 00364 ELSE IF( SIN.GT.STMP ) THEN 00365 VMAX = SIN / STMP 00366 ELSE IF( SIN.LT.EPS*STMP ) THEN 00367 VMAX = ONE / EPS 00368 ELSE IF( SIN.LT.STMP ) THEN 00369 VMAX = STMP / SIN 00370 ELSE 00371 VMAX = ONE 00372 END IF 00373 IF( VMAX.GT.RMAX( 3 ) ) THEN 00374 RMAX( 3 ) = VMAX 00375 IF( NINFO( 3 ).EQ.0 ) 00376 $ LMAX( 3 ) = KNT 00377 END IF 00378 * 00379 * Compare condition numbers for invariant subspace 00380 * without taking its condition number into account 00381 * 00382 IF( SEPIN.LE.V .AND. SEPTMP.LE.V ) THEN 00383 VMAX = ONE 00384 ELSE IF( EPS*SEPIN.GT.SEPTMP ) THEN 00385 VMAX = ONE / EPS 00386 ELSE IF( SEPIN.GT.SEPTMP ) THEN 00387 VMAX = SEPIN / SEPTMP 00388 ELSE IF( SEPIN.LT.EPS*SEPTMP ) THEN 00389 VMAX = ONE / EPS 00390 ELSE IF( SEPIN.LT.SEPTMP ) THEN 00391 VMAX = SEPTMP / SEPIN 00392 ELSE 00393 VMAX = ONE 00394 END IF 00395 IF( VMAX.GT.RMAX( 3 ) ) THEN 00396 RMAX( 3 ) = VMAX 00397 IF( NINFO( 3 ).EQ.0 ) 00398 $ LMAX( 3 ) = KNT 00399 END IF 00400 * 00401 * Compute eigenvalue condition number only and compare 00402 * Update Q 00403 * 00404 VMAX = ZERO 00405 CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00406 CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00407 SEPTMP = -ONE 00408 STMP = -ONE 00409 CALL ZTRSEN( 'E', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00410 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00411 IF( INFO.NE.0 ) THEN 00412 LMAX( 3 ) = KNT 00413 NINFO( 3 ) = NINFO( 3 ) + 1 00414 GO TO 200 00415 END IF 00416 IF( S.NE.STMP ) 00417 $ VMAX = ONE / EPS 00418 IF( -ONE.NE.SEPTMP ) 00419 $ VMAX = ONE / EPS 00420 DO 130 I = 1, N 00421 DO 120 J = 1, N 00422 IF( TTMP( I, J ).NE.T( I, J ) ) 00423 $ VMAX = ONE / EPS 00424 IF( QTMP( I, J ).NE.Q( I, J ) ) 00425 $ VMAX = ONE / EPS 00426 120 CONTINUE 00427 130 CONTINUE 00428 * 00429 * Compute invariant subspace condition number only and compare 00430 * Update Q 00431 * 00432 CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00433 CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00434 SEPTMP = -ONE 00435 STMP = -ONE 00436 CALL ZTRSEN( 'V', 'V', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00437 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00438 IF( INFO.NE.0 ) THEN 00439 LMAX( 3 ) = KNT 00440 NINFO( 3 ) = NINFO( 3 ) + 1 00441 GO TO 200 00442 END IF 00443 IF( -ONE.NE.STMP ) 00444 $ VMAX = ONE / EPS 00445 IF( SEP.NE.SEPTMP ) 00446 $ VMAX = ONE / EPS 00447 DO 150 I = 1, N 00448 DO 140 J = 1, N 00449 IF( TTMP( I, J ).NE.T( I, J ) ) 00450 $ VMAX = ONE / EPS 00451 IF( QTMP( I, J ).NE.Q( I, J ) ) 00452 $ VMAX = ONE / EPS 00453 140 CONTINUE 00454 150 CONTINUE 00455 * 00456 * Compute eigenvalue condition number only and compare 00457 * Do not update Q 00458 * 00459 CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00460 CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00461 SEPTMP = -ONE 00462 STMP = -ONE 00463 CALL ZTRSEN( 'E', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00464 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00465 IF( INFO.NE.0 ) THEN 00466 LMAX( 3 ) = KNT 00467 NINFO( 3 ) = NINFO( 3 ) + 1 00468 GO TO 200 00469 END IF 00470 IF( S.NE.STMP ) 00471 $ VMAX = ONE / EPS 00472 IF( -ONE.NE.SEPTMP ) 00473 $ VMAX = ONE / EPS 00474 DO 170 I = 1, N 00475 DO 160 J = 1, N 00476 IF( TTMP( I, J ).NE.T( I, J ) ) 00477 $ VMAX = ONE / EPS 00478 IF( QTMP( I, J ).NE.QSAV( I, J ) ) 00479 $ VMAX = ONE / EPS 00480 160 CONTINUE 00481 170 CONTINUE 00482 * 00483 * Compute invariant subspace condition number only and compare 00484 * Do not update Q 00485 * 00486 CALL ZLACPY( 'F', N, N, TSAV1, LDT, TTMP, LDT ) 00487 CALL ZLACPY( 'F', N, N, QSAV, LDT, QTMP, LDT ) 00488 SEPTMP = -ONE 00489 STMP = -ONE 00490 CALL ZTRSEN( 'V', 'N', SELECT, N, TTMP, LDT, QTMP, LDT, WTMP, 00491 $ M, STMP, SEPTMP, WORK, LWORK, INFO ) 00492 IF( INFO.NE.0 ) THEN 00493 LMAX( 3 ) = KNT 00494 NINFO( 3 ) = NINFO( 3 ) + 1 00495 GO TO 200 00496 END IF 00497 IF( -ONE.NE.STMP ) 00498 $ VMAX = ONE / EPS 00499 IF( SEP.NE.SEPTMP ) 00500 $ VMAX = ONE / EPS 00501 DO 190 I = 1, N 00502 DO 180 J = 1, N 00503 IF( TTMP( I, J ).NE.T( I, J ) ) 00504 $ VMAX = ONE / EPS 00505 IF( QTMP( I, J ).NE.QSAV( I, J ) ) 00506 $ VMAX = ONE / EPS 00507 180 CONTINUE 00508 190 CONTINUE 00509 IF( VMAX.GT.RMAX( 1 ) ) THEN 00510 RMAX( 1 ) = VMAX 00511 IF( NINFO( 1 ).EQ.0 ) 00512 $ LMAX( 1 ) = KNT 00513 END IF 00514 200 CONTINUE 00515 GO TO 10 00516 * 00517 * End of ZGET38 00518 * 00519 END