LAPACK: TESTING/EIG/dget36.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGET36( 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, LMAX, NIN
00009       DOUBLE PRECISION   RMAX
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            NINFO( 3 )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DGET36 tests DTREXC, a routine for moving blocks (either 1 by 1 or
00019 *  2 by 2) on the diagonal of a matrix in real Schur form.  Thus, DLAEXC
00020 *  computes an orthogonal matrix Q such that
00021 *
00022 *     Q' * T1 * Q  = T2
00023 *
00024 *  and where one of the diagonal blocks of T1 (the one at row IFST) has
00025 *  been moved to position ILST.
00026 *
00027 *  The test code verifies that the residual Q'*T1*Q-T2 is small, that T2
00028 *  is in Schur form, and that the final position of the IFST block is
00029 *  ILST (within +-1).
00030 *
00031 *  The test matrices are read from a file with logical unit number NIN.
00032 *
00033 *  Arguments
00034 *  ==========
00035 *
00036 *  RMAX    (output) DOUBLE PRECISION
00037 *          Value of the largest test ratio.
00038 *
00039 *  LMAX    (output) INTEGER
00040 *          Example number where largest test ratio achieved.
00041 *
00042 *  NINFO   (output) INTEGER array, dimension (3)
00043 *          NINFO(J) is the number of examples where INFO=J.
00044 *
00045 *  KNT     (output) INTEGER
00046 *          Total number of examples tested.
00047 *
00048 *  NIN     (input) INTEGER
00049 *          Input logical unit number.
00050 *
00051 *  =====================================================================
00052 *
00053 *     .. Parameters ..
00054       DOUBLE PRECISION   ZERO, ONE
00055       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00056       INTEGER            LDT, LWORK
00057       PARAMETER          ( LDT = 10, LWORK = 2*LDT*LDT )
00058 *     ..
00059 *     .. Local Scalars ..
00060       INTEGER            I, IFST, IFST1, IFST2, IFSTSV, ILST, ILST1,
00061      $                   ILST2, ILSTSV, INFO1, INFO2, J, LOC, N
00062       DOUBLE PRECISION   EPS, RES
00063 *     ..
00064 *     .. Local Arrays ..
00065       DOUBLE PRECISION   Q( LDT, LDT ), RESULT( 2 ), T1( LDT, LDT ),
00066      $                   T2( LDT, LDT ), TMP( LDT, LDT ), WORK( LWORK )
00067 *     ..
00068 *     .. External Functions ..
00069       DOUBLE PRECISION   DLAMCH
00070       EXTERNAL           DLAMCH
00071 *     ..
00072 *     .. External Subroutines ..
00073       EXTERNAL           DHST01, DLACPY, DLASET, DTREXC
00074 *     ..
00075 *     .. Intrinsic Functions ..
00076       INTRINSIC          ABS, SIGN
00077 *     ..
00078 *     .. Executable Statements ..
00079 *
00080       EPS = DLAMCH( 'P' )
00081       RMAX = ZERO
00082       LMAX = 0
00083       KNT = 0
00084       NINFO( 1 ) = 0
00085       NINFO( 2 ) = 0
00086       NINFO( 3 ) = 0
00087 *
00088 *     Read input data until N=0
00089 *
00090    10 CONTINUE
00091       READ( NIN, FMT = * )N, IFST, ILST
00092       IF( N.EQ.0 )
00093      $   RETURN
00094       KNT = KNT + 1
00095       DO 20 I = 1, N
00096          READ( NIN, FMT = * )( TMP( I, J ), J = 1, N )
00097    20 CONTINUE
00098       CALL DLACPY( 'F', N, N, TMP, LDT, T1, LDT )
00099       CALL DLACPY( 'F', N, N, TMP, LDT, T2, LDT )
00100       IFSTSV = IFST
00101       ILSTSV = ILST
00102       IFST1 = IFST
00103       ILST1 = ILST
00104       IFST2 = IFST
00105       ILST2 = ILST
00106       RES = ZERO
00107 *
00108 *     Test without accumulating Q
00109 *
00110       CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDT )
00111       CALL DTREXC( 'N', N, T1, LDT, Q, LDT, IFST1, ILST1, WORK, INFO1 )
00112       DO 40 I = 1, N
00113          DO 30 J = 1, N
00114             IF( I.EQ.J .AND. Q( I, J ).NE.ONE )
00115      $         RES = RES + ONE / EPS
00116             IF( I.NE.J .AND. Q( I, J ).NE.ZERO )
00117      $         RES = RES + ONE / EPS
00118    30    CONTINUE
00119    40 CONTINUE
00120 *
00121 *     Test with accumulating Q
00122 *
00123       CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDT )
00124       CALL DTREXC( 'V', N, T2, LDT, Q, LDT, IFST2, ILST2, WORK, INFO2 )
00125 *
00126 *     Compare T1 with T2
00127 *
00128       DO 60 I = 1, N
00129          DO 50 J = 1, N
00130             IF( T1( I, J ).NE.T2( I, J ) )
00131      $         RES = RES + ONE / EPS
00132    50    CONTINUE
00133    60 CONTINUE
00134       IF( IFST1.NE.IFST2 )
00135      $   RES = RES + ONE / EPS
00136       IF( ILST1.NE.ILST2 )
00137      $   RES = RES + ONE / EPS
00138       IF( INFO1.NE.INFO2 )
00139      $   RES = RES + ONE / EPS
00140 *
00141 *     Test for successful reordering of T2
00142 *
00143       IF( INFO2.NE.0 ) THEN
00144          NINFO( INFO2 ) = NINFO( INFO2 ) + 1
00145       ELSE
00146          IF( ABS( IFST2-IFSTSV ).GT.1 )
00147      $      RES = RES + ONE / EPS
00148          IF( ABS( ILST2-ILSTSV ).GT.1 )
00149      $      RES = RES + ONE / EPS
00150       END IF
00151 *
00152 *     Test for small residual, and orthogonality of Q
00153 *
00154       CALL DHST01( N, 1, N, TMP, LDT, T2, LDT, Q, LDT, WORK, LWORK,
00155      $             RESULT )
00156       RES = RES + RESULT( 1 ) + RESULT( 2 )
00157 *
00158 *     Test for T2 being in Schur form
00159 *
00160       LOC = 1
00161    70 CONTINUE
00162       IF( T2( LOC+1, LOC ).NE.ZERO ) THEN
00163 *
00164 *        2 by 2 block
00165 *
00166          IF( T2( LOC, LOC+1 ).EQ.ZERO .OR. T2( LOC, LOC ).NE.
00167      $       T2( LOC+1, LOC+1 ) .OR. SIGN( ONE, T2( LOC, LOC+1 ) ).EQ.
00168      $       SIGN( ONE, T2( LOC+1, LOC ) ) )RES = RES + ONE / EPS
00169          DO 80 I = LOC + 2, N
00170             IF( T2( I, LOC ).NE.ZERO )
00171      $         RES = RES + ONE / RES
00172             IF( T2( I, LOC+1 ).NE.ZERO )
00173      $         RES = RES + ONE / RES
00174    80    CONTINUE
00175          LOC = LOC + 2
00176       ELSE
00177 *
00178 *        1 by 1 block
00179 *
00180          DO 90 I = LOC + 1, N
00181             IF( T2( I, LOC ).NE.ZERO )
00182      $         RES = RES + ONE / RES
00183    90    CONTINUE
00184          LOC = LOC + 1
00185       END IF
00186       IF( LOC.LT.N )
00187      $   GO TO 70
00188       IF( RES.GT.RMAX ) THEN
00189          RMAX = RES
00190          LMAX = KNT
00191       END IF
00192       GO TO 10
00193 *
00194 *     End of DGET36
00195 *
00196       END