![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHBT21 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 ZHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, 00012 * RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER UPLO 00016 * INTEGER KA, KS, LDA, LDU, N 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 00020 * COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> ZHBT21 generally checks a decomposition of the form 00030 *> 00031 *> A = U S UC> 00032 *> where * means conjugate transpose, A is hermitian banded, U is 00033 *> unitary, and S is diagonal (if KS=0) or symmetric 00034 *> tridiagonal (if KS=1). 00035 *> 00036 *> Specifically: 00037 *> 00038 *> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp ) 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] UPLO 00045 *> \verbatim 00046 *> UPLO is CHARACTER 00047 *> If UPLO='U', the upper triangle of A and V will be used and 00048 *> the (strictly) lower triangle will not be referenced. 00049 *> If UPLO='L', the lower triangle of A and V will be used and 00050 *> the (strictly) upper triangle will not be referenced. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] N 00054 *> \verbatim 00055 *> N is INTEGER 00056 *> The size of the matrix. If it is zero, ZHBT21 does nothing. 00057 *> It must be at least zero. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] KA 00061 *> \verbatim 00062 *> KA is INTEGER 00063 *> The bandwidth of the matrix A. It must be at least zero. If 00064 *> it is larger than N-1, then max( 0, N-1 ) will be used. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] KS 00068 *> \verbatim 00069 *> KS is INTEGER 00070 *> The bandwidth of the matrix S. It may only be zero or one. 00071 *> If zero, then S is diagonal, and E is not referenced. If 00072 *> one, then S is symmetric tri-diagonal. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] A 00076 *> \verbatim 00077 *> A is COMPLEX*16 array, dimension (LDA, N) 00078 *> The original (unfactored) matrix. It is assumed to be 00079 *> hermitian, and only the upper (UPLO='U') or only the lower 00080 *> (UPLO='L') will be referenced. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] LDA 00084 *> \verbatim 00085 *> LDA is INTEGER 00086 *> The leading dimension of A. It must be at least 1 00087 *> and at least min( KA, N-1 ). 00088 *> \endverbatim 00089 *> 00090 *> \param[in] D 00091 *> \verbatim 00092 *> D is DOUBLE PRECISION array, dimension (N) 00093 *> The diagonal of the (symmetric tri-) diagonal matrix S. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] E 00097 *> \verbatim 00098 *> E is DOUBLE PRECISION array, dimension (N-1) 00099 *> The off-diagonal of the (symmetric tri-) diagonal matrix S. 00100 *> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 00101 *> (3,2) element, etc. 00102 *> Not referenced if KS=0. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] U 00106 *> \verbatim 00107 *> U is COMPLEX*16 array, dimension (LDU, N) 00108 *> The unitary matrix in the decomposition, expressed as a 00109 *> dense matrix (i.e., not as a product of Householder 00110 *> transformations, Givens transformations, etc.) 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LDU 00114 *> \verbatim 00115 *> LDU is INTEGER 00116 *> The leading dimension of U. LDU must be at least N and 00117 *> at least 1. 00118 *> \endverbatim 00119 *> 00120 *> \param[out] WORK 00121 *> \verbatim 00122 *> WORK is COMPLEX*16 array, dimension (N**2) 00123 *> \endverbatim 00124 *> 00125 *> \param[out] RWORK 00126 *> \verbatim 00127 *> RWORK is DOUBLE PRECISION array, dimension (N) 00128 *> \endverbatim 00129 *> 00130 *> \param[out] RESULT 00131 *> \verbatim 00132 *> RESULT is DOUBLE PRECISION array, dimension (2) 00133 *> The values computed by the two tests described above. The 00134 *> values are currently limited to 1/ulp, to avoid overflow. 00135 *> \endverbatim 00136 * 00137 * Authors: 00138 * ======== 00139 * 00140 *> \author Univ. of Tennessee 00141 *> \author Univ. of California Berkeley 00142 *> \author Univ. of Colorado Denver 00143 *> \author NAG Ltd. 00144 * 00145 *> \date November 2011 00146 * 00147 *> \ingroup complex16_eig 00148 * 00149 * ===================================================================== 00150 SUBROUTINE ZHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, 00151 $ RWORK, RESULT ) 00152 * 00153 * -- LAPACK test routine (version 3.4.0) -- 00154 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00155 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00156 * November 2011 00157 * 00158 * .. Scalar Arguments .. 00159 CHARACTER UPLO 00160 INTEGER KA, KS, LDA, LDU, N 00161 * .. 00162 * .. Array Arguments .. 00163 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 00164 COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * ) 00165 * .. 00166 * 00167 * ===================================================================== 00168 * 00169 * .. Parameters .. 00170 COMPLEX*16 CZERO, CONE 00171 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00172 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00173 DOUBLE PRECISION ZERO, ONE 00174 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00175 * .. 00176 * .. Local Scalars .. 00177 LOGICAL LOWER 00178 CHARACTER CUPLO 00179 INTEGER IKA, J, JC, JR 00180 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM 00181 * .. 00182 * .. External Functions .. 00183 LOGICAL LSAME 00184 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHB, ZLANHP 00185 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHB, ZLANHP 00186 * .. 00187 * .. External Subroutines .. 00188 EXTERNAL ZGEMM, ZHPR, ZHPR2 00189 * .. 00190 * .. Intrinsic Functions .. 00191 INTRINSIC DBLE, DCMPLX, MAX, MIN 00192 * .. 00193 * .. Executable Statements .. 00194 * 00195 * Constants 00196 * 00197 RESULT( 1 ) = ZERO 00198 RESULT( 2 ) = ZERO 00199 IF( N.LE.0 ) 00200 $ RETURN 00201 * 00202 IKA = MAX( 0, MIN( N-1, KA ) ) 00203 * 00204 IF( LSAME( UPLO, 'U' ) ) THEN 00205 LOWER = .FALSE. 00206 CUPLO = 'U' 00207 ELSE 00208 LOWER = .TRUE. 00209 CUPLO = 'L' 00210 END IF 00211 * 00212 UNFL = DLAMCH( 'Safe minimum' ) 00213 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00214 * 00215 * Some Error Checks 00216 * 00217 * Do Test 1 00218 * 00219 * Norm of A: 00220 * 00221 ANORM = MAX( ZLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL ) 00222 * 00223 * Compute error matrix: Error = A - U S U* 00224 * 00225 * Copy A from SB to SP storage format. 00226 * 00227 J = 0 00228 DO 50 JC = 1, N 00229 IF( LOWER ) THEN 00230 DO 10 JR = 1, MIN( IKA+1, N+1-JC ) 00231 J = J + 1 00232 WORK( J ) = A( JR, JC ) 00233 10 CONTINUE 00234 DO 20 JR = IKA + 2, N + 1 - JC 00235 J = J + 1 00236 WORK( J ) = ZERO 00237 20 CONTINUE 00238 ELSE 00239 DO 30 JR = IKA + 2, JC 00240 J = J + 1 00241 WORK( J ) = ZERO 00242 30 CONTINUE 00243 DO 40 JR = MIN( IKA, JC-1 ), 0, -1 00244 J = J + 1 00245 WORK( J ) = A( IKA+1-JR, JC ) 00246 40 CONTINUE 00247 END IF 00248 50 CONTINUE 00249 * 00250 DO 60 J = 1, N 00251 CALL ZHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK ) 00252 60 CONTINUE 00253 * 00254 IF( N.GT.1 .AND. KS.EQ.1 ) THEN 00255 DO 70 J = 1, N - 1 00256 CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1, 00257 $ U( 1, J+1 ), 1, WORK ) 00258 70 CONTINUE 00259 END IF 00260 WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK ) 00261 * 00262 IF( ANORM.GT.WNORM ) THEN 00263 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 00264 ELSE 00265 IF( ANORM.LT.ONE ) THEN 00266 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 00267 ELSE 00268 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 00269 END IF 00270 END IF 00271 * 00272 * Do Test 2 00273 * 00274 * Compute UU* - I 00275 * 00276 CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK, 00277 $ N ) 00278 * 00279 DO 80 J = 1, N 00280 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE 00281 80 CONTINUE 00282 * 00283 RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ), 00284 $ DBLE( N ) ) / ( N*ULP ) 00285 * 00286 RETURN 00287 * 00288 * End of ZHBT21 00289 * 00290 END