1*> \brief \b ZHBT21 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, 12* RWORK, RESULT ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER KA, KS, LDA, LDU, N 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 20* COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> ZHBT21 generally checks a decomposition of the form 30*> 31*> A = U S U**H 32*> 33*> where **H means conjugate transpose, A is hermitian banded, U is 34*> unitary, and S is diagonal (if KS=0) or symmetric 35*> tridiagonal (if KS=1). 36*> 37*> Specifically: 38*> 39*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) and 40*> RESULT(2) = | I - U U**H | / ( n ulp ) 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] UPLO 47*> \verbatim 48*> UPLO is CHARACTER 49*> If UPLO='U', the upper triangle of A and V will be used and 50*> the (strictly) lower triangle will not be referenced. 51*> If UPLO='L', the lower triangle of A and V will be used and 52*> the (strictly) upper triangle will not be referenced. 53*> \endverbatim 54*> 55*> \param[in] N 56*> \verbatim 57*> N is INTEGER 58*> The size of the matrix. If it is zero, ZHBT21 does nothing. 59*> It must be at least zero. 60*> \endverbatim 61*> 62*> \param[in] KA 63*> \verbatim 64*> KA is INTEGER 65*> The bandwidth of the matrix A. It must be at least zero. If 66*> it is larger than N-1, then max( 0, N-1 ) will be used. 67*> \endverbatim 68*> 69*> \param[in] KS 70*> \verbatim 71*> KS is INTEGER 72*> The bandwidth of the matrix S. It may only be zero or one. 73*> If zero, then S is diagonal, and E is not referenced. If 74*> one, then S is symmetric tri-diagonal. 75*> \endverbatim 76*> 77*> \param[in] A 78*> \verbatim 79*> A is COMPLEX*16 array, dimension (LDA, N) 80*> The original (unfactored) matrix. It is assumed to be 81*> hermitian, and only the upper (UPLO='U') or only the lower 82*> (UPLO='L') will be referenced. 83*> \endverbatim 84*> 85*> \param[in] LDA 86*> \verbatim 87*> LDA is INTEGER 88*> The leading dimension of A. It must be at least 1 89*> and at least min( KA, N-1 ). 90*> \endverbatim 91*> 92*> \param[in] D 93*> \verbatim 94*> D is DOUBLE PRECISION array, dimension (N) 95*> The diagonal of the (symmetric tri-) diagonal matrix S. 96*> \endverbatim 97*> 98*> \param[in] E 99*> \verbatim 100*> E is DOUBLE PRECISION array, dimension (N-1) 101*> The off-diagonal of the (symmetric tri-) diagonal matrix S. 102*> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 103*> (3,2) element, etc. 104*> Not referenced if KS=0. 105*> \endverbatim 106*> 107*> \param[in] U 108*> \verbatim 109*> U is COMPLEX*16 array, dimension (LDU, N) 110*> The unitary matrix in the decomposition, expressed as a 111*> dense matrix (i.e., not as a product of Householder 112*> transformations, Givens transformations, etc.) 113*> \endverbatim 114*> 115*> \param[in] LDU 116*> \verbatim 117*> LDU is INTEGER 118*> The leading dimension of U. LDU must be at least N and 119*> at least 1. 120*> \endverbatim 121*> 122*> \param[out] WORK 123*> \verbatim 124*> WORK is COMPLEX*16 array, dimension (N**2) 125*> \endverbatim 126*> 127*> \param[out] RWORK 128*> \verbatim 129*> RWORK is DOUBLE PRECISION array, dimension (N) 130*> \endverbatim 131*> 132*> \param[out] RESULT 133*> \verbatim 134*> RESULT is DOUBLE PRECISION array, dimension (2) 135*> The values computed by the two tests described above. The 136*> values are currently limited to 1/ulp, to avoid overflow. 137*> \endverbatim 138* 139* Authors: 140* ======== 141* 142*> \author Univ. of Tennessee 143*> \author Univ. of California Berkeley 144*> \author Univ. of Colorado Denver 145*> \author NAG Ltd. 146* 147*> \ingroup complex16_eig 148* 149* ===================================================================== 150 SUBROUTINE ZHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, 151 $ RWORK, RESULT ) 152* 153* -- LAPACK test routine -- 154* -- LAPACK is a software package provided by Univ. of Tennessee, -- 155* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 156* 157* .. Scalar Arguments .. 158 CHARACTER UPLO 159 INTEGER KA, KS, LDA, LDU, N 160* .. 161* .. Array Arguments .. 162 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 163 COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * ) 164* .. 165* 166* ===================================================================== 167* 168* .. Parameters .. 169 COMPLEX*16 CZERO, CONE 170 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 171 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 172 DOUBLE PRECISION ZERO, ONE 173 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 174* .. 175* .. Local Scalars .. 176 LOGICAL LOWER 177 CHARACTER CUPLO 178 INTEGER IKA, J, JC, JR 179 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM 180* .. 181* .. External Functions .. 182 LOGICAL LSAME 183 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHB, ZLANHP 184 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHB, ZLANHP 185* .. 186* .. External Subroutines .. 187 EXTERNAL ZGEMM, ZHPR, ZHPR2 188* .. 189* .. Intrinsic Functions .. 190 INTRINSIC DBLE, DCMPLX, MAX, MIN 191* .. 192* .. Executable Statements .. 193* 194* Constants 195* 196 RESULT( 1 ) = ZERO 197 RESULT( 2 ) = ZERO 198 IF( N.LE.0 ) 199 $ RETURN 200* 201 IKA = MAX( 0, MIN( N-1, KA ) ) 202* 203 IF( LSAME( UPLO, 'U' ) ) THEN 204 LOWER = .FALSE. 205 CUPLO = 'U' 206 ELSE 207 LOWER = .TRUE. 208 CUPLO = 'L' 209 END IF 210* 211 UNFL = DLAMCH( 'Safe minimum' ) 212 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 213* 214* Some Error Checks 215* 216* Do Test 1 217* 218* Norm of A: 219* 220 ANORM = MAX( ZLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL ) 221* 222* Compute error matrix: Error = A - U S U**H 223* 224* Copy A from SB to SP storage format. 225* 226 J = 0 227 DO 50 JC = 1, N 228 IF( LOWER ) THEN 229 DO 10 JR = 1, MIN( IKA+1, N+1-JC ) 230 J = J + 1 231 WORK( J ) = A( JR, JC ) 232 10 CONTINUE 233 DO 20 JR = IKA + 2, N + 1 - JC 234 J = J + 1 235 WORK( J ) = ZERO 236 20 CONTINUE 237 ELSE 238 DO 30 JR = IKA + 2, JC 239 J = J + 1 240 WORK( J ) = ZERO 241 30 CONTINUE 242 DO 40 JR = MIN( IKA, JC-1 ), 0, -1 243 J = J + 1 244 WORK( J ) = A( IKA+1-JR, JC ) 245 40 CONTINUE 246 END IF 247 50 CONTINUE 248* 249 DO 60 J = 1, N 250 CALL ZHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK ) 251 60 CONTINUE 252* 253 IF( N.GT.1 .AND. KS.EQ.1 ) THEN 254 DO 70 J = 1, N - 1 255 CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1, 256 $ U( 1, J+1 ), 1, WORK ) 257 70 CONTINUE 258 END IF 259 WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK ) 260* 261 IF( ANORM.GT.WNORM ) THEN 262 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 263 ELSE 264 IF( ANORM.LT.ONE ) THEN 265 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 266 ELSE 267 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 268 END IF 269 END IF 270* 271* Do Test 2 272* 273* Compute U U**H - I 274* 275 CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK, 276 $ N ) 277* 278 DO 80 J = 1, N 279 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE 280 80 CONTINUE 281* 282 RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ), 283 $ DBLE( N ) ) / ( N*ULP ) 284* 285 RETURN 286* 287* End of ZHBT21 288* 289 END 290