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