LAPACK: TESTING/EIG/cget38.f Source File

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