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