![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGET36 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 DGET36( RMAX, LMAX, NINFO, KNT, NIN ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER KNT, LMAX, NIN 00015 * DOUBLE PRECISION RMAX 00016 * .. 00017 * .. Array Arguments .. 00018 * INTEGER NINFO( 3 ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> DGET36 tests DTREXC, a routine for moving blocks (either 1 by 1 or 00028 *> 2 by 2) on the diagonal of a matrix in real Schur form. Thus, DLAEXC 00029 *> computes an orthogonal matrix Q such that 00030 *> 00031 *> Q' * T1 * Q = T2 00032 *> 00033 *> and where one of the diagonal blocks of T1 (the one at row IFST) has 00034 *> been moved to position ILST. 00035 *> 00036 *> The test code verifies that the residual Q'*T1*Q-T2 is small, that T2 00037 *> is in Schur form, and that the final position of the IFST block is 00038 *> ILST (within +-1). 00039 *> 00040 *> The test matrices are read from a file with logical unit number NIN. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[out] RMAX 00047 *> \verbatim 00048 *> RMAX is DOUBLE PRECISION 00049 *> Value of the largest test ratio. 00050 *> \endverbatim 00051 *> 00052 *> \param[out] LMAX 00053 *> \verbatim 00054 *> LMAX is INTEGER 00055 *> Example number where largest test ratio achieved. 00056 *> \endverbatim 00057 *> 00058 *> \param[out] NINFO 00059 *> \verbatim 00060 *> NINFO is INTEGER array, dimension (3) 00061 *> NINFO(J) is the number of examples where INFO=J. 00062 *> \endverbatim 00063 *> 00064 *> \param[out] KNT 00065 *> \verbatim 00066 *> KNT is INTEGER 00067 *> Total number of examples tested. 00068 *> \endverbatim 00069 *> 00070 *> \param[in] NIN 00071 *> \verbatim 00072 *> NIN is INTEGER 00073 *> Input logical unit number. 00074 *> \endverbatim 00075 * 00076 * Authors: 00077 * ======== 00078 * 00079 *> \author Univ. of Tennessee 00080 *> \author Univ. of California Berkeley 00081 *> \author Univ. of Colorado Denver 00082 *> \author NAG Ltd. 00083 * 00084 *> \date November 2011 00085 * 00086 *> \ingroup double_eig 00087 * 00088 * ===================================================================== 00089 SUBROUTINE DGET36( RMAX, LMAX, NINFO, KNT, NIN ) 00090 * 00091 * -- LAPACK test routine (version 3.4.0) -- 00092 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00093 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00094 * November 2011 00095 * 00096 * .. Scalar Arguments .. 00097 INTEGER KNT, LMAX, NIN 00098 DOUBLE PRECISION RMAX 00099 * .. 00100 * .. Array Arguments .. 00101 INTEGER NINFO( 3 ) 00102 * .. 00103 * 00104 * ===================================================================== 00105 * 00106 * .. Parameters .. 00107 DOUBLE PRECISION ZERO, ONE 00108 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00109 INTEGER LDT, LWORK 00110 PARAMETER ( LDT = 10, LWORK = 2*LDT*LDT ) 00111 * .. 00112 * .. Local Scalars .. 00113 INTEGER I, IFST, IFST1, IFST2, IFSTSV, ILST, ILST1, 00114 $ ILST2, ILSTSV, INFO1, INFO2, J, LOC, N 00115 DOUBLE PRECISION EPS, RES 00116 * .. 00117 * .. Local Arrays .. 00118 DOUBLE PRECISION Q( LDT, LDT ), RESULT( 2 ), T1( LDT, LDT ), 00119 $ T2( LDT, LDT ), TMP( LDT, LDT ), WORK( LWORK ) 00120 * .. 00121 * .. External Functions .. 00122 DOUBLE PRECISION DLAMCH 00123 EXTERNAL DLAMCH 00124 * .. 00125 * .. External Subroutines .. 00126 EXTERNAL DHST01, DLACPY, DLASET, DTREXC 00127 * .. 00128 * .. Intrinsic Functions .. 00129 INTRINSIC ABS, SIGN 00130 * .. 00131 * .. Executable Statements .. 00132 * 00133 EPS = DLAMCH( 'P' ) 00134 RMAX = ZERO 00135 LMAX = 0 00136 KNT = 0 00137 NINFO( 1 ) = 0 00138 NINFO( 2 ) = 0 00139 NINFO( 3 ) = 0 00140 * 00141 * Read input data until N=0 00142 * 00143 10 CONTINUE 00144 READ( NIN, FMT = * )N, IFST, ILST 00145 IF( N.EQ.0 ) 00146 $ RETURN 00147 KNT = KNT + 1 00148 DO 20 I = 1, N 00149 READ( NIN, FMT = * )( TMP( I, J ), J = 1, N ) 00150 20 CONTINUE 00151 CALL DLACPY( 'F', N, N, TMP, LDT, T1, LDT ) 00152 CALL DLACPY( 'F', N, N, TMP, LDT, T2, LDT ) 00153 IFSTSV = IFST 00154 ILSTSV = ILST 00155 IFST1 = IFST 00156 ILST1 = ILST 00157 IFST2 = IFST 00158 ILST2 = ILST 00159 RES = ZERO 00160 * 00161 * Test without accumulating Q 00162 * 00163 CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDT ) 00164 CALL DTREXC( 'N', N, T1, LDT, Q, LDT, IFST1, ILST1, WORK, INFO1 ) 00165 DO 40 I = 1, N 00166 DO 30 J = 1, N 00167 IF( I.EQ.J .AND. Q( I, J ).NE.ONE ) 00168 $ RES = RES + ONE / EPS 00169 IF( I.NE.J .AND. Q( I, J ).NE.ZERO ) 00170 $ RES = RES + ONE / EPS 00171 30 CONTINUE 00172 40 CONTINUE 00173 * 00174 * Test with accumulating Q 00175 * 00176 CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDT ) 00177 CALL DTREXC( 'V', N, T2, LDT, Q, LDT, IFST2, ILST2, WORK, INFO2 ) 00178 * 00179 * Compare T1 with T2 00180 * 00181 DO 60 I = 1, N 00182 DO 50 J = 1, N 00183 IF( T1( I, J ).NE.T2( I, J ) ) 00184 $ RES = RES + ONE / EPS 00185 50 CONTINUE 00186 60 CONTINUE 00187 IF( IFST1.NE.IFST2 ) 00188 $ RES = RES + ONE / EPS 00189 IF( ILST1.NE.ILST2 ) 00190 $ RES = RES + ONE / EPS 00191 IF( INFO1.NE.INFO2 ) 00192 $ RES = RES + ONE / EPS 00193 * 00194 * Test for successful reordering of T2 00195 * 00196 IF( INFO2.NE.0 ) THEN 00197 NINFO( INFO2 ) = NINFO( INFO2 ) + 1 00198 ELSE 00199 IF( ABS( IFST2-IFSTSV ).GT.1 ) 00200 $ RES = RES + ONE / EPS 00201 IF( ABS( ILST2-ILSTSV ).GT.1 ) 00202 $ RES = RES + ONE / EPS 00203 END IF 00204 * 00205 * Test for small residual, and orthogonality of Q 00206 * 00207 CALL DHST01( N, 1, N, TMP, LDT, T2, LDT, Q, LDT, WORK, LWORK, 00208 $ RESULT ) 00209 RES = RES + RESULT( 1 ) + RESULT( 2 ) 00210 * 00211 * Test for T2 being in Schur form 00212 * 00213 LOC = 1 00214 70 CONTINUE 00215 IF( T2( LOC+1, LOC ).NE.ZERO ) THEN 00216 * 00217 * 2 by 2 block 00218 * 00219 IF( T2( LOC, LOC+1 ).EQ.ZERO .OR. T2( LOC, LOC ).NE. 00220 $ T2( LOC+1, LOC+1 ) .OR. SIGN( ONE, T2( LOC, LOC+1 ) ).EQ. 00221 $ SIGN( ONE, T2( LOC+1, LOC ) ) )RES = RES + ONE / EPS 00222 DO 80 I = LOC + 2, N 00223 IF( T2( I, LOC ).NE.ZERO ) 00224 $ RES = RES + ONE / RES 00225 IF( T2( I, LOC+1 ).NE.ZERO ) 00226 $ RES = RES + ONE / RES 00227 80 CONTINUE 00228 LOC = LOC + 2 00229 ELSE 00230 * 00231 * 1 by 1 block 00232 * 00233 DO 90 I = LOC + 1, N 00234 IF( T2( I, LOC ).NE.ZERO ) 00235 $ RES = RES + ONE / RES 00236 90 CONTINUE 00237 LOC = LOC + 1 00238 END IF 00239 IF( LOC.LT.N ) 00240 $ GO TO 70 00241 IF( RES.GT.RMAX ) THEN 00242 RMAX = RES 00243 LMAX = KNT 00244 END IF 00245 GO TO 10 00246 * 00247 * End of DGET36 00248 * 00249 END