LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zget38.f
Go to the documentation of this file.
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
 All Files Functions