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