1*> \brief \b ZHETRI2X 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHETRI2X + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri2x.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri2x.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri2x.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, LDA, N, NB 26* .. 27* .. Array Arguments .. 28* INTEGER IPIV( * ) 29* COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix 39*> A using the factorization A = U*D*U**H or A = L*D*L**H computed by 40*> ZHETRF. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] UPLO 47*> \verbatim 48*> UPLO is CHARACTER*1 49*> Specifies whether the details of the factorization are stored 50*> as an upper or lower triangular matrix. 51*> = 'U': Upper triangular, form is A = U*D*U**H; 52*> = 'L': Lower triangular, form is A = L*D*L**H. 53*> \endverbatim 54*> 55*> \param[in] N 56*> \verbatim 57*> N is INTEGER 58*> The order of the matrix A. N >= 0. 59*> \endverbatim 60*> 61*> \param[in,out] A 62*> \verbatim 63*> A is COMPLEX*16 array, dimension (LDA,N) 64*> On entry, the NNB diagonal matrix D and the multipliers 65*> used to obtain the factor U or L as computed by ZHETRF. 66*> 67*> On exit, if INFO = 0, the (symmetric) inverse of the original 68*> matrix. If UPLO = 'U', the upper triangular part of the 69*> inverse is formed and the part of A below the diagonal is not 70*> referenced; if UPLO = 'L' the lower triangular part of the 71*> inverse is formed and the part of A above the diagonal is 72*> not referenced. 73*> \endverbatim 74*> 75*> \param[in] LDA 76*> \verbatim 77*> LDA is INTEGER 78*> The leading dimension of the array A. LDA >= max(1,N). 79*> \endverbatim 80*> 81*> \param[in] IPIV 82*> \verbatim 83*> IPIV is INTEGER array, dimension (N) 84*> Details of the interchanges and the NNB structure of D 85*> as determined by ZHETRF. 86*> \endverbatim 87*> 88*> \param[out] WORK 89*> \verbatim 90*> WORK is COMPLEX*16 array, dimension (N+NB+1,NB+3) 91*> \endverbatim 92*> 93*> \param[in] NB 94*> \verbatim 95*> NB is INTEGER 96*> Block size 97*> \endverbatim 98*> 99*> \param[out] INFO 100*> \verbatim 101*> INFO is INTEGER 102*> = 0: successful exit 103*> < 0: if INFO = -i, the i-th argument had an illegal value 104*> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its 105*> inverse could not be computed. 106*> \endverbatim 107* 108* Authors: 109* ======== 110* 111*> \author Univ. of Tennessee 112*> \author Univ. of California Berkeley 113*> \author Univ. of Colorado Denver 114*> \author NAG Ltd. 115* 116*> \ingroup complex16HEcomputational 117* 118* ===================================================================== 119 SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) 120* 121* -- LAPACK computational routine -- 122* -- LAPACK is a software package provided by Univ. of Tennessee, -- 123* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 124* 125* .. Scalar Arguments .. 126 CHARACTER UPLO 127 INTEGER INFO, LDA, N, NB 128* .. 129* .. Array Arguments .. 130 INTEGER IPIV( * ) 131 COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) 132* .. 133* 134* ===================================================================== 135* 136* .. Parameters .. 137 DOUBLE PRECISION ONE 138 COMPLEX*16 CONE, ZERO 139 PARAMETER ( ONE = 1.0D+0, 140 $ CONE = ( 1.0D+0, 0.0D+0 ), 141 $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 142* .. 143* .. Local Scalars .. 144 LOGICAL UPPER 145 INTEGER I, IINFO, IP, K, CUT, NNB 146 INTEGER COUNT 147 INTEGER J, U11, INVD 148 149 COMPLEX*16 AK, AKKP1, AKP1, D, T 150 COMPLEX*16 U01_I_J, U01_IP1_J 151 COMPLEX*16 U11_I_J, U11_IP1_J 152* .. 153* .. External Functions .. 154 LOGICAL LSAME 155 EXTERNAL LSAME 156* .. 157* .. External Subroutines .. 158 EXTERNAL ZSYCONV, XERBLA, ZTRTRI 159 EXTERNAL ZGEMM, ZTRMM, ZHESWAPR 160* .. 161* .. Intrinsic Functions .. 162 INTRINSIC MAX 163* .. 164* .. Executable Statements .. 165* 166* Test the input parameters. 167* 168 INFO = 0 169 UPPER = LSAME( UPLO, 'U' ) 170 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 171 INFO = -1 172 ELSE IF( N.LT.0 ) THEN 173 INFO = -2 174 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 175 INFO = -4 176 END IF 177* 178* Quick return if possible 179* 180* 181 IF( INFO.NE.0 ) THEN 182 CALL XERBLA( 'ZHETRI2X', -INFO ) 183 RETURN 184 END IF 185 IF( N.EQ.0 ) 186 $ RETURN 187* 188* Convert A 189* Workspace got Non-diag elements of D 190* 191 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 192* 193* Check that the diagonal matrix D is nonsingular. 194* 195 IF( UPPER ) THEN 196* 197* Upper triangular storage: examine D from bottom to top 198* 199 DO INFO = N, 1, -1 200 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 201 $ RETURN 202 END DO 203 ELSE 204* 205* Lower triangular storage: examine D from top to bottom. 206* 207 DO INFO = 1, N 208 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 209 $ RETURN 210 END DO 211 END IF 212 INFO = 0 213* 214* Splitting Workspace 215* U01 is a block (N,NB+1) 216* The first element of U01 is in WORK(1,1) 217* U11 is a block (NB+1,NB+1) 218* The first element of U11 is in WORK(N+1,1) 219 U11 = N 220* INVD is a block (N,2) 221* The first element of INVD is in WORK(1,INVD) 222 INVD = NB+2 223 224 IF( UPPER ) THEN 225* 226* invA = P * inv(U**H)*inv(D)*inv(U)*P**H. 227* 228 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) 229* 230* inv(D) and inv(D)*inv(U) 231* 232 K=1 233 DO WHILE ( K .LE. N ) 234 IF( IPIV( K ).GT.0 ) THEN 235* 1 x 1 diagonal NNB 236 WORK(K,INVD) = ONE / REAL ( A( K, K ) ) 237 WORK(K,INVD+1) = 0 238 K=K+1 239 ELSE 240* 2 x 2 diagonal NNB 241 T = ABS ( WORK(K+1,1) ) 242 AK = DBLE ( A( K, K ) ) / T 243 AKP1 = DBLE ( A( K+1, K+1 ) ) / T 244 AKKP1 = WORK(K+1,1) / T 245 D = T*( AK*AKP1-ONE ) 246 WORK(K,INVD) = AKP1 / D 247 WORK(K+1,INVD+1) = AK / D 248 WORK(K,INVD+1) = -AKKP1 / D 249 WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) ) 250 K=K+2 251 END IF 252 END DO 253* 254* inv(U**H) = (inv(U))**H 255* 256* inv(U**H)*inv(D)*inv(U) 257* 258 CUT=N 259 DO WHILE (CUT .GT. 0) 260 NNB=NB 261 IF (CUT .LE. NNB) THEN 262 NNB=CUT 263 ELSE 264 COUNT = 0 265* count negative elements, 266 DO I=CUT+1-NNB,CUT 267 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 268 END DO 269* need a even number for a clear cut 270 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 271 END IF 272 273 CUT=CUT-NNB 274* 275* U01 Block 276* 277 DO I=1,CUT 278 DO J=1,NNB 279 WORK(I,J)=A(I,CUT+J) 280 END DO 281 END DO 282* 283* U11 Block 284* 285 DO I=1,NNB 286 WORK(U11+I,I)=CONE 287 DO J=1,I-1 288 WORK(U11+I,J)=ZERO 289 END DO 290 DO J=I+1,NNB 291 WORK(U11+I,J)=A(CUT+I,CUT+J) 292 END DO 293 END DO 294* 295* invD*U01 296* 297 I=1 298 DO WHILE (I .LE. CUT) 299 IF (IPIV(I) > 0) THEN 300 DO J=1,NNB 301 WORK(I,J)=WORK(I,INVD)*WORK(I,J) 302 END DO 303 I=I+1 304 ELSE 305 DO J=1,NNB 306 U01_I_J = WORK(I,J) 307 U01_IP1_J = WORK(I+1,J) 308 WORK(I,J)=WORK(I,INVD)*U01_I_J+ 309 $ WORK(I,INVD+1)*U01_IP1_J 310 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+ 311 $ WORK(I+1,INVD+1)*U01_IP1_J 312 END DO 313 I=I+2 314 END IF 315 END DO 316* 317* invD1*U11 318* 319 I=1 320 DO WHILE (I .LE. NNB) 321 IF (IPIV(CUT+I) > 0) THEN 322 DO J=I,NNB 323 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 324 END DO 325 I=I+1 326 ELSE 327 DO J=I,NNB 328 U11_I_J = WORK(U11+I,J) 329 U11_IP1_J = WORK(U11+I+1,J) 330 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 331 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J) 332 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+ 333 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J 334 END DO 335 I=I+2 336 END IF 337 END DO 338* 339* U11**H*invD1*U11->U11 340* 341 CALL ZTRMM('L','U','C','U',NNB, NNB, 342 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 343* 344 DO I=1,NNB 345 DO J=I,NNB 346 A(CUT+I,CUT+J)=WORK(U11+I,J) 347 END DO 348 END DO 349* 350* U01**H*invD*U01->A(CUT+I,CUT+J) 351* 352 CALL ZGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA, 353 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 354* 355* U11 = U11**H*invD1*U11 + U01**H*invD*U01 356* 357 DO I=1,NNB 358 DO J=I,NNB 359 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 360 END DO 361 END DO 362* 363* U01 = U00**H*invD0*U01 364* 365 CALL ZTRMM('L',UPLO,'C','U',CUT, NNB, 366 $ CONE,A,LDA,WORK,N+NB+1) 367 368* 369* Update U01 370* 371 DO I=1,CUT 372 DO J=1,NNB 373 A(I,CUT+J)=WORK(I,J) 374 END DO 375 END DO 376* 377* Next Block 378* 379 END DO 380* 381* Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H 382* 383 I=1 384 DO WHILE ( I .LE. N ) 385 IF( IPIV(I) .GT. 0 ) THEN 386 IP=IPIV(I) 387 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP ) 388 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I ) 389 ELSE 390 IP=-IPIV(I) 391 I=I+1 392 IF ( (I-1) .LT. IP) 393 $ CALL ZHESWAPR( UPLO, N, A, LDA, I-1 ,IP ) 394 IF ( (I-1) .GT. IP) 395 $ CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I-1 ) 396 ENDIF 397 I=I+1 398 END DO 399 ELSE 400* 401* LOWER... 402* 403* invA = P * inv(U**H)*inv(D)*inv(U)*P**H. 404* 405 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) 406* 407* inv(D) and inv(D)*inv(U) 408* 409 K=N 410 DO WHILE ( K .GE. 1 ) 411 IF( IPIV( K ).GT.0 ) THEN 412* 1 x 1 diagonal NNB 413 WORK(K,INVD) = ONE / REAL ( A( K, K ) ) 414 WORK(K,INVD+1) = 0 415 K=K-1 416 ELSE 417* 2 x 2 diagonal NNB 418 T = ABS ( WORK(K-1,1) ) 419 AK = DBLE ( A( K-1, K-1 ) ) / T 420 AKP1 = DBLE ( A( K, K ) ) / T 421 AKKP1 = WORK(K-1,1) / T 422 D = T*( AK*AKP1-ONE ) 423 WORK(K-1,INVD) = AKP1 / D 424 WORK(K,INVD) = AK / D 425 WORK(K,INVD+1) = -AKKP1 / D 426 WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) ) 427 K=K-2 428 END IF 429 END DO 430* 431* inv(U**H) = (inv(U))**H 432* 433* inv(U**H)*inv(D)*inv(U) 434* 435 CUT=0 436 DO WHILE (CUT .LT. N) 437 NNB=NB 438 IF (CUT + NNB .GE. N) THEN 439 NNB=N-CUT 440 ELSE 441 COUNT = 0 442* count negative elements, 443 DO I=CUT+1,CUT+NNB 444 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 445 END DO 446* need a even number for a clear cut 447 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 448 END IF 449* L21 Block 450 DO I=1,N-CUT-NNB 451 DO J=1,NNB 452 WORK(I,J)=A(CUT+NNB+I,CUT+J) 453 END DO 454 END DO 455* L11 Block 456 DO I=1,NNB 457 WORK(U11+I,I)=CONE 458 DO J=I+1,NNB 459 WORK(U11+I,J)=ZERO 460 END DO 461 DO J=1,I-1 462 WORK(U11+I,J)=A(CUT+I,CUT+J) 463 END DO 464 END DO 465* 466* invD*L21 467* 468 I=N-CUT-NNB 469 DO WHILE (I .GE. 1) 470 IF (IPIV(CUT+NNB+I) > 0) THEN 471 DO J=1,NNB 472 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 473 END DO 474 I=I-1 475 ELSE 476 DO J=1,NNB 477 U01_I_J = WORK(I,J) 478 U01_IP1_J = WORK(I-1,J) 479 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 480 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 481 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 482 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 483 END DO 484 I=I-2 485 END IF 486 END DO 487* 488* invD1*L11 489* 490 I=NNB 491 DO WHILE (I .GE. 1) 492 IF (IPIV(CUT+I) > 0) THEN 493 DO J=1,NNB 494 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 495 END DO 496 I=I-1 497 ELSE 498 DO J=1,NNB 499 U11_I_J = WORK(U11+I,J) 500 U11_IP1_J = WORK(U11+I-1,J) 501 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 502 $ WORK(CUT+I,INVD+1)*U11_IP1_J 503 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 504 $ WORK(CUT+I-1,INVD)*U11_IP1_J 505 END DO 506 I=I-2 507 END IF 508 END DO 509* 510* L11**H*invD1*L11->L11 511* 512 CALL ZTRMM('L',UPLO,'C','U',NNB, NNB, 513 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 514* 515 DO I=1,NNB 516 DO J=1,I 517 A(CUT+I,CUT+J)=WORK(U11+I,J) 518 END DO 519 END DO 520* 521 IF ( (CUT+NNB) .LT. N ) THEN 522* 523* L21**H*invD2*L21->A(CUT+I,CUT+J) 524* 525 CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1) 526 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 527 528* 529* L11 = L11**H*invD1*L11 + U01**H*invD*U01 530* 531 DO I=1,NNB 532 DO J=1,I 533 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 534 END DO 535 END DO 536* 537* L01 = L22**H*invD2*L21 538* 539 CALL ZTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB, 540 $ CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 541 542* Update L21 543 DO I=1,N-CUT-NNB 544 DO J=1,NNB 545 A(CUT+NNB+I,CUT+J)=WORK(I,J) 546 END DO 547 END DO 548 ELSE 549* 550* L11 = L11**H*invD1*L11 551* 552 DO I=1,NNB 553 DO J=1,I 554 A(CUT+I,CUT+J)=WORK(U11+I,J) 555 END DO 556 END DO 557 END IF 558* 559* Next Block 560* 561 CUT=CUT+NNB 562 END DO 563* 564* Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H 565* 566 I=N 567 DO WHILE ( I .GE. 1 ) 568 IF( IPIV(I) .GT. 0 ) THEN 569 IP=IPIV(I) 570 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP ) 571 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I ) 572 ELSE 573 IP=-IPIV(I) 574 IF ( I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP ) 575 IF ( I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I ) 576 I=I-1 577 ENDIF 578 I=I-1 579 END DO 580 END IF 581* 582 RETURN 583* 584* End of ZHETRI2X 585* 586 END 587 588