LAPACK: TESTING/EIG/dget38.f Source File

Go to the documentation of this file.
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