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