![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSTT22 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 DSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, 00012 * LDWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER KBAND, LDU, LDWORK, M, N 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), SD( * ), 00019 * $ SE( * ), U( LDU, * ), WORK( LDWORK, * ) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> DSTT22 checks a set of M eigenvalues and eigenvectors, 00029 *> 00030 *> A U = U S 00031 *> 00032 *> where A is symmetric tridiagonal, the columns of U are orthogonal, 00033 *> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1). 00034 *> Two tests are performed: 00035 *> 00036 *> RESULT(1) = | U' A U - S | / ( |A| m ulp ) 00037 *> 00038 *> RESULT(2) = | I - U'U | / ( m ulp ) 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] N 00045 *> \verbatim 00046 *> N is INTEGER 00047 *> The size of the matrix. If it is zero, DSTT22 does nothing. 00048 *> It must be at least zero. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] M 00052 *> \verbatim 00053 *> M is INTEGER 00054 *> The number of eigenpairs to check. If it is zero, DSTT22 00055 *> does nothing. It must be at least zero. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] KBAND 00059 *> \verbatim 00060 *> KBAND is INTEGER 00061 *> The bandwidth of the matrix S. It may only be zero or one. 00062 *> If zero, then S is diagonal, and SE is not referenced. If 00063 *> one, then S is symmetric tri-diagonal. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] AD 00067 *> \verbatim 00068 *> AD is DOUBLE PRECISION array, dimension (N) 00069 *> The diagonal of the original (unfactored) matrix A. A is 00070 *> assumed to be symmetric tridiagonal. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] AE 00074 *> \verbatim 00075 *> AE is DOUBLE PRECISION array, dimension (N) 00076 *> The off-diagonal of the original (unfactored) matrix A. A 00077 *> is assumed to be symmetric tridiagonal. AE(1) is ignored, 00078 *> AE(2) is the (1,2) and (2,1) element, etc. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] SD 00082 *> \verbatim 00083 *> SD is DOUBLE PRECISION array, dimension (N) 00084 *> The diagonal of the (symmetric tri-) diagonal matrix S. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] SE 00088 *> \verbatim 00089 *> SE is DOUBLE PRECISION array, dimension (N) 00090 *> The off-diagonal of the (symmetric tri-) diagonal matrix S. 00091 *> Not referenced if KBSND=0. If KBAND=1, then AE(1) is 00092 *> ignored, SE(2) is the (1,2) and (2,1) element, etc. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] U 00096 *> \verbatim 00097 *> U is DOUBLE PRECISION array, dimension (LDU, N) 00098 *> The orthogonal matrix in the decomposition. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDU 00102 *> \verbatim 00103 *> LDU is INTEGER 00104 *> The leading dimension of U. LDU must be at least N. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] WORK 00108 *> \verbatim 00109 *> WORK is DOUBLE PRECISION array, dimension (LDWORK, M+1) 00110 *> \endverbatim 00111 *> 00112 *> \param[in] LDWORK 00113 *> \verbatim 00114 *> LDWORK is INTEGER 00115 *> The leading dimension of WORK. LDWORK must be at least 00116 *> max(1,M). 00117 *> \endverbatim 00118 *> 00119 *> \param[out] RESULT 00120 *> \verbatim 00121 *> RESULT is DOUBLE PRECISION array, dimension (2) 00122 *> The values computed by the two tests described above. The 00123 *> values are currently limited to 1/ulp, to avoid overflow. 00124 *> \endverbatim 00125 * 00126 * Authors: 00127 * ======== 00128 * 00129 *> \author Univ. of Tennessee 00130 *> \author Univ. of California Berkeley 00131 *> \author Univ. of Colorado Denver 00132 *> \author NAG Ltd. 00133 * 00134 *> \date November 2011 00135 * 00136 *> \ingroup double_eig 00137 * 00138 * ===================================================================== 00139 SUBROUTINE DSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, 00140 $ LDWORK, RESULT ) 00141 * 00142 * -- LAPACK test routine (version 3.4.0) -- 00143 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00144 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00145 * November 2011 00146 * 00147 * .. Scalar Arguments .. 00148 INTEGER KBAND, LDU, LDWORK, M, N 00149 * .. 00150 * .. Array Arguments .. 00151 DOUBLE PRECISION AD( * ), AE( * ), RESULT( 2 ), SD( * ), 00152 $ SE( * ), U( LDU, * ), WORK( LDWORK, * ) 00153 * .. 00154 * 00155 * ===================================================================== 00156 * 00157 * .. Parameters .. 00158 DOUBLE PRECISION ZERO, ONE 00159 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00160 * .. 00161 * .. Local Scalars .. 00162 INTEGER I, J, K 00163 DOUBLE PRECISION ANORM, AUKJ, ULP, UNFL, WNORM 00164 * .. 00165 * .. External Functions .. 00166 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00167 EXTERNAL DLAMCH, DLANGE, DLANSY 00168 * .. 00169 * .. External Subroutines .. 00170 EXTERNAL DGEMM 00171 * .. 00172 * .. Intrinsic Functions .. 00173 INTRINSIC ABS, DBLE, MAX, MIN 00174 * .. 00175 * .. Executable Statements .. 00176 * 00177 RESULT( 1 ) = ZERO 00178 RESULT( 2 ) = ZERO 00179 IF( N.LE.0 .OR. M.LE.0 ) 00180 $ RETURN 00181 * 00182 UNFL = DLAMCH( 'Safe minimum' ) 00183 ULP = DLAMCH( 'Epsilon' ) 00184 * 00185 * Do Test 1 00186 * 00187 * Compute the 1-norm of A. 00188 * 00189 IF( N.GT.1 ) THEN 00190 ANORM = ABS( AD( 1 ) ) + ABS( AE( 1 ) ) 00191 DO 10 J = 2, N - 1 00192 ANORM = MAX( ANORM, ABS( AD( J ) )+ABS( AE( J ) )+ 00193 $ ABS( AE( J-1 ) ) ) 00194 10 CONTINUE 00195 ANORM = MAX( ANORM, ABS( AD( N ) )+ABS( AE( N-1 ) ) ) 00196 ELSE 00197 ANORM = ABS( AD( 1 ) ) 00198 END IF 00199 ANORM = MAX( ANORM, UNFL ) 00200 * 00201 * Norm of U'AU - S 00202 * 00203 DO 40 I = 1, M 00204 DO 30 J = 1, M 00205 WORK( I, J ) = ZERO 00206 DO 20 K = 1, N 00207 AUKJ = AD( K )*U( K, J ) 00208 IF( K.NE.N ) 00209 $ AUKJ = AUKJ + AE( K )*U( K+1, J ) 00210 IF( K.NE.1 ) 00211 $ AUKJ = AUKJ + AE( K-1 )*U( K-1, J ) 00212 WORK( I, J ) = WORK( I, J ) + U( K, I )*AUKJ 00213 20 CONTINUE 00214 30 CONTINUE 00215 WORK( I, I ) = WORK( I, I ) - SD( I ) 00216 IF( KBAND.EQ.1 ) THEN 00217 IF( I.NE.1 ) 00218 $ WORK( I, I-1 ) = WORK( I, I-1 ) - SE( I-1 ) 00219 IF( I.NE.N ) 00220 $ WORK( I, I+1 ) = WORK( I, I+1 ) - SE( I ) 00221 END IF 00222 40 CONTINUE 00223 * 00224 WNORM = DLANSY( '1', 'L', M, WORK, M, WORK( 1, M+1 ) ) 00225 * 00226 IF( ANORM.GT.WNORM ) THEN 00227 RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP ) 00228 ELSE 00229 IF( ANORM.LT.ONE ) THEN 00230 RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP ) 00231 ELSE 00232 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( M ) ) / ( M*ULP ) 00233 END IF 00234 END IF 00235 * 00236 * Do Test 2 00237 * 00238 * Compute U'U - I 00239 * 00240 CALL DGEMM( 'T', 'N', M, M, N, ONE, U, LDU, U, LDU, ZERO, WORK, 00241 $ M ) 00242 * 00243 DO 50 J = 1, M 00244 WORK( J, J ) = WORK( J, J ) - ONE 00245 50 CONTINUE 00246 * 00247 RESULT( 2 ) = MIN( DBLE( M ), DLANGE( '1', M, M, WORK, M, WORK( 1, 00248 $ M+1 ) ) ) / ( M*ULP ) 00249 * 00250 RETURN 00251 * 00252 * End of DSTT22 00253 * 00254 END