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