1*> \brief \b ZHPT21 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 ZHPT21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, 12* TAU, WORK, RWORK, RESULT ) 13* 14* .. Scalar Arguments .. 15* CHARACTER UPLO 16* INTEGER ITYPE, KBAND, LDU, N 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 20* COMPLEX*16 AP( * ), TAU( * ), U( LDU, * ), VP( * ), 21* $ WORK( * ) 22* .. 23* 24* 25*> \par Purpose: 26* ============= 27*> 28*> \verbatim 29*> 30*> ZHPT21 generally checks a decomposition of the form 31*> 32*> A = U S UC> 33*> where * means conjugate transpose, A is hermitian, U is 34*> unitary, and S is diagonal (if KBAND=0) or (real) symmetric 35*> tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as 36*> a dense matrix, otherwise the U is expressed as a product of 37*> Householder transformations, whose vectors are stored in the 38*> array "V" and whose scaling constants are in "TAU"; we shall 39*> use the letter "V" to refer to the product of Householder 40*> transformations (which should be equal to U). 41*> 42*> Specifically, if ITYPE=1, then: 43*> 44*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp ) 45*> 46*> If ITYPE=2, then: 47*> 48*> RESULT(1) = | A - V S V* | / ( |A| n ulp ) 49*> 50*> If ITYPE=3, then: 51*> 52*> RESULT(1) = | I - UV* | / ( n ulp ) 53*> 54*> Packed storage means that, for example, if UPLO='U', then the columns 55*> of the upper triangle of A are stored one after another, so that 56*> A(1,j+1) immediately follows A(j,j) in the array AP. Similarly, if 57*> UPLO='L', then the columns of the lower triangle of A are stored one 58*> after another in AP, so that A(j+1,j+1) immediately follows A(n,j) 59*> in the array AP. This means that A(i,j) is stored in: 60*> 61*> AP( i + j*(j-1)/2 ) if UPLO='U' 62*> 63*> AP( i + (2*n-j)*(j-1)/2 ) if UPLO='L' 64*> 65*> The array VP bears the same relation to the matrix V that A does to 66*> AP. 67*> 68*> For ITYPE > 1, the transformation U is expressed as a product 69*> of Householder transformations: 70*> 71*> If UPLO='U', then V = H(n-1)...H(1), where 72*> 73*> H(j) = I - tau(j) v(j) v(j)C> 74*> and the first j-1 elements of v(j) are stored in V(1:j-1,j+1), 75*> (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ), 76*> the j-th element is 1, and the last n-j elements are 0. 77*> 78*> If UPLO='L', then V = H(1)...H(n-1), where 79*> 80*> H(j) = I - tau(j) v(j) v(j)C> 81*> and the first j elements of v(j) are 0, the (j+1)-st is 1, and the 82*> (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e., 83*> in VP( (2*n-j)*(j-1)/2 + j+2 : (2*n-j)*(j-1)/2 + n ) .) 84*> \endverbatim 85* 86* Arguments: 87* ========== 88* 89*> \param[in] ITYPE 90*> \verbatim 91*> ITYPE is INTEGER 92*> Specifies the type of tests to be performed. 93*> 1: U expressed as a dense unitary matrix: 94*> RESULT(1) = | A - U S U* | / ( |A| n ulp ) *andC> RESULT(2) = | I - UU* | / ( n ulp ) 95*> 96*> 2: U expressed as a product V of Housholder transformations: 97*> RESULT(1) = | A - V S V* | / ( |A| n ulp ) 98*> 99*> 3: U expressed both as a dense unitary matrix and 100*> as a product of Housholder transformations: 101*> RESULT(1) = | I - UV* | / ( n ulp ) 102*> \endverbatim 103*> 104*> \param[in] UPLO 105*> \verbatim 106*> UPLO is CHARACTER 107*> If UPLO='U', the upper triangle of A and V will be used and 108*> the (strictly) lower triangle will not be referenced. 109*> If UPLO='L', the lower triangle of A and V will be used and 110*> the (strictly) upper triangle will not be referenced. 111*> \endverbatim 112*> 113*> \param[in] N 114*> \verbatim 115*> N is INTEGER 116*> The size of the matrix. If it is zero, ZHPT21 does nothing. 117*> It must be at least zero. 118*> \endverbatim 119*> 120*> \param[in] KBAND 121*> \verbatim 122*> KBAND is INTEGER 123*> The bandwidth of the matrix. It may only be zero or one. 124*> If zero, then S is diagonal, and E is not referenced. If 125*> one, then S is symmetric tri-diagonal. 126*> \endverbatim 127*> 128*> \param[in] AP 129*> \verbatim 130*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 131*> The original (unfactored) matrix. It is assumed to be 132*> hermitian, and contains the columns of just the upper 133*> triangle (UPLO='U') or only the lower triangle (UPLO='L'), 134*> packed one after another. 135*> \endverbatim 136*> 137*> \param[in] D 138*> \verbatim 139*> D is DOUBLE PRECISION array, dimension (N) 140*> The diagonal of the (symmetric tri-) diagonal matrix. 141*> \endverbatim 142*> 143*> \param[in] E 144*> \verbatim 145*> E is DOUBLE PRECISION array, dimension (N) 146*> The off-diagonal of the (symmetric tri-) diagonal matrix. 147*> E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and 148*> (3,2) element, etc. 149*> Not referenced if KBAND=0. 150*> \endverbatim 151*> 152*> \param[in] U 153*> \verbatim 154*> U is COMPLEX*16 array, dimension (LDU, N) 155*> If ITYPE=1 or 3, this contains the unitary matrix in 156*> the decomposition, expressed as a dense matrix. If ITYPE=2, 157*> then it is not referenced. 158*> \endverbatim 159*> 160*> \param[in] LDU 161*> \verbatim 162*> LDU is INTEGER 163*> The leading dimension of U. LDU must be at least N and 164*> at least 1. 165*> \endverbatim 166*> 167*> \param[in] VP 168*> \verbatim 169*> VP is DOUBLE PRECISION array, dimension (N*(N+1)/2) 170*> If ITYPE=2 or 3, the columns of this array contain the 171*> Householder vectors used to describe the unitary matrix 172*> in the decomposition, as described in purpose. 173*> *NOTE* If ITYPE=2 or 3, V is modified and restored. The 174*> subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U') 175*> is set to one, and later reset to its original value, during 176*> the course of the calculation. 177*> If ITYPE=1, then it is neither referenced nor modified. 178*> \endverbatim 179*> 180*> \param[in] TAU 181*> \verbatim 182*> TAU is COMPLEX*16 array, dimension (N) 183*> If ITYPE >= 2, then TAU(j) is the scalar factor of 184*> v(j) v(j)* in the Householder transformation H(j) of 185*> the product U = H(1)...H(n-2) 186*> If ITYPE < 2, then TAU is not referenced. 187*> \endverbatim 188*> 189*> \param[out] WORK 190*> \verbatim 191*> WORK is COMPLEX*16 array, dimension (N**2) 192*> Workspace. 193*> \endverbatim 194*> 195*> \param[out] RWORK 196*> \verbatim 197*> RWORK is DOUBLE PRECISION array, dimension (N) 198*> Workspace. 199*> \endverbatim 200*> 201*> \param[out] RESULT 202*> \verbatim 203*> RESULT is DOUBLE PRECISION array, dimension (2) 204*> The values computed by the two tests described above. The 205*> values are currently limited to 1/ulp, to avoid overflow. 206*> RESULT(1) is always modified. RESULT(2) is modified only 207*> if ITYPE=1. 208*> \endverbatim 209* 210* Authors: 211* ======== 212* 213*> \author Univ. of Tennessee 214*> \author Univ. of California Berkeley 215*> \author Univ. of Colorado Denver 216*> \author NAG Ltd. 217* 218*> \date November 2011 219* 220*> \ingroup complex16_eig 221* 222* ===================================================================== 223 SUBROUTINE ZHPT21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, 224 $ TAU, WORK, RWORK, RESULT ) 225* 226* -- LAPACK test routine (version 3.4.0) -- 227* -- LAPACK is a software package provided by Univ. of Tennessee, -- 228* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 229* November 2011 230* 231* .. Scalar Arguments .. 232 CHARACTER UPLO 233 INTEGER ITYPE, KBAND, LDU, N 234* .. 235* .. Array Arguments .. 236 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * ) 237 COMPLEX*16 AP( * ), TAU( * ), U( LDU, * ), VP( * ), 238 $ WORK( * ) 239* .. 240* 241* ===================================================================== 242* 243* .. Parameters .. 244 DOUBLE PRECISION ZERO, ONE, TEN 245 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 10.0D+0 ) 246 DOUBLE PRECISION HALF 247 PARAMETER ( HALF = 1.0D+0 / 2.0D+0 ) 248 COMPLEX*16 CZERO, CONE 249 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 250 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 251* .. 252* .. Local Scalars .. 253 LOGICAL LOWER 254 CHARACTER CUPLO 255 INTEGER IINFO, J, JP, JP1, JR, LAP 256 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM 257 COMPLEX*16 TEMP, VSAVE 258* .. 259* .. External Functions .. 260 LOGICAL LSAME 261 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHP 262 COMPLEX*16 ZDOTC 263 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHP, ZDOTC 264* .. 265* .. External Subroutines .. 266 EXTERNAL ZAXPY, ZCOPY, ZGEMM, ZHPMV, ZHPR, ZHPR2, 267 $ ZLACPY, ZLASET, ZUPMTR 268* .. 269* .. Intrinsic Functions .. 270 INTRINSIC DBLE, DCMPLX, MAX, MIN 271* .. 272* .. Executable Statements .. 273* 274* Constants 275* 276 RESULT( 1 ) = ZERO 277 IF( ITYPE.EQ.1 ) 278 $ RESULT( 2 ) = ZERO 279 IF( N.LE.0 ) 280 $ RETURN 281* 282 LAP = ( N*( N+1 ) ) / 2 283* 284 IF( LSAME( UPLO, 'U' ) ) THEN 285 LOWER = .FALSE. 286 CUPLO = 'U' 287 ELSE 288 LOWER = .TRUE. 289 CUPLO = 'L' 290 END IF 291* 292 UNFL = DLAMCH( 'Safe minimum' ) 293 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 294* 295* Some Error Checks 296* 297 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 298 RESULT( 1 ) = TEN / ULP 299 RETURN 300 END IF 301* 302* Do Test 1 303* 304* Norm of A: 305* 306 IF( ITYPE.EQ.3 ) THEN 307 ANORM = ONE 308 ELSE 309 ANORM = MAX( ZLANHP( '1', CUPLO, N, AP, RWORK ), UNFL ) 310 END IF 311* 312* Compute error matrix: 313* 314 IF( ITYPE.EQ.1 ) THEN 315* 316* ITYPE=1: error = A - U S U* 317* 318 CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N ) 319 CALL ZCOPY( LAP, AP, 1, WORK, 1 ) 320* 321 DO 10 J = 1, N 322 CALL ZHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK ) 323 10 CONTINUE 324* 325 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN 326 DO 20 J = 1, N - 1 327 CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1, 328 $ U( 1, J-1 ), 1, WORK ) 329 20 CONTINUE 330 END IF 331 WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK ) 332* 333 ELSE IF( ITYPE.EQ.2 ) THEN 334* 335* ITYPE=2: error = V S V* - A 336* 337 CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N ) 338* 339 IF( LOWER ) THEN 340 WORK( LAP ) = D( N ) 341 DO 40 J = N - 1, 1, -1 342 JP = ( ( 2*N-J )*( J-1 ) ) / 2 343 JP1 = JP + N - J 344 IF( KBAND.EQ.1 ) THEN 345 WORK( JP+J+1 ) = ( CONE-TAU( J ) )*E( J ) 346 DO 30 JR = J + 2, N 347 WORK( JP+JR ) = -TAU( J )*E( J )*VP( JP+JR ) 348 30 CONTINUE 349 END IF 350* 351 IF( TAU( J ).NE.CZERO ) THEN 352 VSAVE = VP( JP+J+1 ) 353 VP( JP+J+1 ) = CONE 354 CALL ZHPMV( 'L', N-J, CONE, WORK( JP1+J+1 ), 355 $ VP( JP+J+1 ), 1, CZERO, WORK( LAP+1 ), 1 ) 356 TEMP = -HALF*TAU( J )*ZDOTC( N-J, WORK( LAP+1 ), 1, 357 $ VP( JP+J+1 ), 1 ) 358 CALL ZAXPY( N-J, TEMP, VP( JP+J+1 ), 1, WORK( LAP+1 ), 359 $ 1 ) 360 CALL ZHPR2( 'L', N-J, -TAU( J ), VP( JP+J+1 ), 1, 361 $ WORK( LAP+1 ), 1, WORK( JP1+J+1 ) ) 362* 363 VP( JP+J+1 ) = VSAVE 364 END IF 365 WORK( JP+J ) = D( J ) 366 40 CONTINUE 367 ELSE 368 WORK( 1 ) = D( 1 ) 369 DO 60 J = 1, N - 1 370 JP = ( J*( J-1 ) ) / 2 371 JP1 = JP + J 372 IF( KBAND.EQ.1 ) THEN 373 WORK( JP1+J ) = ( CONE-TAU( J ) )*E( J ) 374 DO 50 JR = 1, J - 1 375 WORK( JP1+JR ) = -TAU( J )*E( J )*VP( JP1+JR ) 376 50 CONTINUE 377 END IF 378* 379 IF( TAU( J ).NE.CZERO ) THEN 380 VSAVE = VP( JP1+J ) 381 VP( JP1+J ) = CONE 382 CALL ZHPMV( 'U', J, CONE, WORK, VP( JP1+1 ), 1, CZERO, 383 $ WORK( LAP+1 ), 1 ) 384 TEMP = -HALF*TAU( J )*ZDOTC( J, WORK( LAP+1 ), 1, 385 $ VP( JP1+1 ), 1 ) 386 CALL ZAXPY( J, TEMP, VP( JP1+1 ), 1, WORK( LAP+1 ), 387 $ 1 ) 388 CALL ZHPR2( 'U', J, -TAU( J ), VP( JP1+1 ), 1, 389 $ WORK( LAP+1 ), 1, WORK ) 390 VP( JP1+J ) = VSAVE 391 END IF 392 WORK( JP1+J+1 ) = D( J+1 ) 393 60 CONTINUE 394 END IF 395* 396 DO 70 J = 1, LAP 397 WORK( J ) = WORK( J ) - AP( J ) 398 70 CONTINUE 399 WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK ) 400* 401 ELSE IF( ITYPE.EQ.3 ) THEN 402* 403* ITYPE=3: error = U V* - I 404* 405 IF( N.LT.2 ) 406 $ RETURN 407 CALL ZLACPY( ' ', N, N, U, LDU, WORK, N ) 408 CALL ZUPMTR( 'R', CUPLO, 'C', N, N, VP, TAU, WORK, N, 409 $ WORK( N**2+1 ), IINFO ) 410 IF( IINFO.NE.0 ) THEN 411 RESULT( 1 ) = TEN / ULP 412 RETURN 413 END IF 414* 415 DO 80 J = 1, N 416 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE 417 80 CONTINUE 418* 419 WNORM = ZLANGE( '1', N, N, WORK, N, RWORK ) 420 END IF 421* 422 IF( ANORM.GT.WNORM ) THEN 423 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP ) 424 ELSE 425 IF( ANORM.LT.ONE ) THEN 426 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP ) 427 ELSE 428 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP ) 429 END IF 430 END IF 431* 432* Do Test 2 433* 434* Compute UU* - I 435* 436 IF( ITYPE.EQ.1 ) THEN 437 CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, 438 $ WORK, N ) 439* 440 DO 90 J = 1, N 441 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE 442 90 CONTINUE 443* 444 RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ), 445 $ DBLE( N ) ) / ( N*ULP ) 446 END IF 447* 448 RETURN 449* 450* End of ZHPT21 451* 452 END 453