1*> \brief \b ZHET22 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 ZHET22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU, 12* V, LDV, TAU, WORK, RWORK, RESULT ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 20* COMPLEX*16 A( LDA, * ), TAU( * ), U( LDU, * ), 21* $ V( LDV, * ), WORK( * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> ZHET22 generally checks a decomposition of the form 31*> 32*> A U = U S 33*> 34*> where A is complex Hermitian, the columns of U are orthonormal, 35*> and S is diagonal (if KBAND=0) or symmetric tridiagonal (if 36*> KBAND=1). If ITYPE=1, then U is represented as a dense matrix, 37*> otherwise the U is expressed as a product of Householder 38*> transformations, whose vectors are stored in the array "V" and 39*> whose scaling constants are in "TAU"; we shall use the letter 40*> "V" to refer to the product of Householder transformations 41*> (which should be equal to U). 42*> 43*> Specifically, if ITYPE=1, then: 44*> 45*> RESULT(1) = | U**H A U - S | / ( |A| m ulp ) and 46*> RESULT(2) = | I - U**H U | / ( m ulp ) 47*> \endverbatim 48* 49* Arguments: 50* ========== 51* 52*> \verbatim 53*> ITYPE INTEGER 54*> Specifies the type of tests to be performed. 55*> 1: U expressed as a dense orthogonal matrix: 56*> RESULT(1) = | A - U S U**H | / ( |A| n ulp ) *and 57*> RESULT(2) = | I - U U**H | / ( n ulp ) 58*> 59*> UPLO CHARACTER 60*> If UPLO='U', the upper triangle of A will be used and the 61*> (strictly) lower triangle will not be referenced. If 62*> UPLO='L', the lower triangle of A will be used and the 63*> (strictly) upper triangle will not be referenced. 64*> Not modified. 65*> 66*> N INTEGER 67*> The size of the matrix. If it is zero, ZHET22 does nothing. 68*> It must be at least zero. 69*> Not modified. 70*> 71*> M INTEGER 72*> The number of columns of U. If it is zero, ZHET22 does 73*> nothing. It must be at least zero. 74*> Not modified. 75*> 76*> KBAND INTEGER 77*> The bandwidth of the matrix. It may only be zero or one. 78*> If zero, then S is diagonal, and E is not referenced. If 79*> one, then S is symmetric tri-diagonal. 80*> Not modified. 81*> 82*> A COMPLEX*16 array, dimension (LDA , N) 83*> The original (unfactored) matrix. It is assumed to be 84*> symmetric, and only the upper (UPLO='U') or only the lower 85*> (UPLO='L') will be referenced. 86*> Not modified. 87*> 88*> LDA INTEGER 89*> The leading dimension of A. It must be at least 1 90*> and at least N. 91*> Not modified. 92*> 93*> D DOUBLE PRECISION array, dimension (N) 94*> The diagonal of the (symmetric tri-) diagonal matrix. 95*> Not modified. 96*> 97*> E DOUBLE PRECISION array, dimension (N) 98*> The off-diagonal of the (symmetric tri-) diagonal matrix. 99*> E(1) is ignored, E(2) is the (1,2) and (2,1) element, etc. 100*> Not referenced if KBAND=0. 101*> Not modified. 102*> 103*> U COMPLEX*16 array, dimension (LDU, N) 104*> If ITYPE=1, this contains the orthogonal matrix in 105*> the decomposition, expressed as a dense matrix. 106*> Not modified. 107*> 108*> LDU INTEGER 109*> The leading dimension of U. LDU must be at least N and 110*> at least 1. 111*> Not modified. 112*> 113*> V COMPLEX*16 array, dimension (LDV, N) 114*> If ITYPE=2 or 3, the lower triangle of this array contains 115*> the Householder vectors used to describe the orthogonal 116*> matrix in the decomposition. If ITYPE=1, then it is not 117*> referenced. 118*> Not modified. 119*> 120*> LDV INTEGER 121*> The leading dimension of V. LDV must be at least N and 122*> at least 1. 123*> Not modified. 124*> 125*> TAU COMPLEX*16 array, dimension (N) 126*> If ITYPE >= 2, then TAU(j) is the scalar factor of 127*> v(j) v(j)**H in the Householder transformation H(j) of 128*> the product U = H(1)...H(n-2) 129*> If ITYPE < 2, then TAU is not referenced. 130*> Not modified. 131*> 132*> WORK COMPLEX*16 array, dimension (2*N**2) 133*> Workspace. 134*> Modified. 135*> 136*> RWORK DOUBLE PRECISION array, dimension (N) 137*> Workspace. 138*> Modified. 139*> 140*> RESULT DOUBLE PRECISION array, dimension (2) 141*> The values computed by the two tests described above. The 142*> values are currently limited to 1/ulp, to avoid overflow. 143*> RESULT(1) is always modified. RESULT(2) is modified only 144*> if LDU is at least N. 145*> Modified. 146*> \endverbatim 147* 148* Authors: 149* ======== 150* 151*> \author Univ. of Tennessee 152*> \author Univ. of California Berkeley 153*> \author Univ. of Colorado Denver 154*> \author NAG Ltd. 155* 156*> \ingroup complex16_eig 157* 158* ===================================================================== 159 SUBROUTINE ZHET22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU, 160 $ V, LDV, TAU, WORK, RWORK, RESULT ) 161* 162* -- LAPACK test routine -- 163* -- LAPACK is a software package provided by Univ. of Tennessee, -- 164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 165* 166* .. Scalar Arguments .. 167 CHARACTER UPLO 168 INTEGER ITYPE, KBAND, LDA, LDU, LDV, M, N 169* .. 170* .. Array Arguments .. 171 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 172 COMPLEX*16 A( LDA, * ), TAU( * ), U( LDU, * ), 173 $ V( LDV, * ), WORK( * ) 174* .. 175* 176* ===================================================================== 177* 178* .. Parameters .. 179 DOUBLE PRECISION ZERO, ONE 180 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 181 COMPLEX*16 CZERO, CONE 182 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), 183 $ CONE = ( 1.0D0, 0.0D0 ) ) 184* .. 185* .. Local Scalars .. 186 INTEGER J, JJ, JJ1, JJ2, NN, NNP1 187 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM 188* .. 189* .. External Functions .. 190 DOUBLE PRECISION DLAMCH, ZLANHE 191 EXTERNAL DLAMCH, ZLANHE 192* .. 193* .. External Subroutines .. 194 EXTERNAL ZGEMM, ZHEMM, ZUNT01 195* .. 196* .. Intrinsic Functions .. 197 INTRINSIC DBLE, MAX, MIN 198* .. 199* .. Executable Statements .. 200* 201 RESULT( 1 ) = ZERO 202 RESULT( 2 ) = ZERO 203 IF( N.LE.0 .OR. M.LE.0 ) 204 $ RETURN 205* 206 UNFL = DLAMCH( 'Safe minimum' ) 207 ULP = DLAMCH( 'Precision' ) 208* 209* Do Test 1 210* 211* Norm of A: 212* 213 ANORM = MAX( ZLANHE( '1', UPLO, N, A, LDA, RWORK ), UNFL ) 214* 215* Compute error matrix: 216* 217* ITYPE=1: error = U**H A U - S 218* 219 CALL ZHEMM( 'L', UPLO, N, M, CONE, A, LDA, U, LDU, CZERO, WORK, 220 $ N ) 221 NN = N*N 222 NNP1 = NN + 1 223 CALL ZGEMM( 'C', 'N', M, M, N, CONE, U, LDU, WORK, N, CZERO, 224 $ WORK( NNP1 ), N ) 225 DO 10 J = 1, M 226 JJ = NN + ( J-1 )*N + J 227 WORK( JJ ) = WORK( JJ ) - D( J ) 228 10 CONTINUE 229 IF( KBAND.EQ.1 .AND. N.GT.1 ) THEN 230 DO 20 J = 2, M 231 JJ1 = NN + ( J-1 )*N + J - 1 232 JJ2 = NN + ( J-2 )*N + J 233 WORK( JJ1 ) = WORK( JJ1 ) - E( J-1 ) 234 WORK( JJ2 ) = WORK( JJ2 ) - E( J-1 ) 235 20 CONTINUE 236 END IF 237 WNORM = ZLANHE( '1', UPLO, M, WORK( NNP1 ), N, RWORK ) 238* 239 IF( ANORM.GT.WNORM ) THEN 240 RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP ) 241 ELSE 242 IF( ANORM.LT.ONE ) THEN 243 RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP ) 244 ELSE 245 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( M ) ) / ( M*ULP ) 246 END IF 247 END IF 248* 249* Do Test 2 250* 251* Compute U**H U - I 252* 253 IF( ITYPE.EQ.1 ) 254 $ CALL ZUNT01( 'Columns', N, M, U, LDU, WORK, 2*N*N, RWORK, 255 $ RESULT( 2 ) ) 256* 257 RETURN 258* 259* End of ZHET22 260* 261 END 262