1*> \brief \b DSYTRI2X 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DSYTRI2X + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri2x.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri2x.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri2x.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSYTRI2X( 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* DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DSYTRI2X computes the inverse of a real symmetric indefinite matrix 39*> A using the factorization A = U*D*U**T or A = L*D*L**T computed by 40*> DSYTRF. 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**T; 52*> = 'L': Lower triangular, form is A = L*D*L**T. 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 DOUBLE PRECISION 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 DSYTRF. 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 DSYTRF. 86*> \endverbatim 87*> 88*> \param[out] WORK 89*> \verbatim 90*> WORK is DOUBLE PRECISION 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 June 2017 117* 118*> \ingroup doubleSYcomputational 119* 120* ===================================================================== 121 SUBROUTINE DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) 122* 123* -- LAPACK computational routine (version 3.7.1) -- 124* -- LAPACK is a software package provided by Univ. of Tennessee, -- 125* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 126* June 2017 127* 128* .. Scalar Arguments .. 129 CHARACTER UPLO 130 INTEGER INFO, LDA, N, NB 131* .. 132* .. Array Arguments .. 133 INTEGER IPIV( * ) 134 DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* ) 135* .. 136* 137* ===================================================================== 138* 139* .. Parameters .. 140 DOUBLE PRECISION ONE, ZERO 141 PARAMETER ( ONE = 1.0D+0, ZERO = 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 DOUBLE PRECISION AK, AKKP1, AKP1, D, T 150 DOUBLE PRECISION U01_I_J, U01_IP1_J 151 DOUBLE PRECISION U11_I_J, U11_IP1_J 152* .. 153* .. External Functions .. 154 LOGICAL LSAME 155 EXTERNAL LSAME 156* .. 157* .. External Subroutines .. 158 EXTERNAL DSYCONV, XERBLA, DTRTRI 159 EXTERNAL DGEMM, DTRMM, DSYSWAPR 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( 'DSYTRI2X', -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 DSYCONV( 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**T)*inv(D)*inv(U)*P**T. 227* 228 CALL DTRTRI( 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 / A( K, K ) 237 WORK(K,INVD+1) = 0 238 K=K+1 239 ELSE 240* 2 x 2 diagonal NNB 241 T = WORK(K+1,1) 242 AK = A( K, K ) / T 243 AKP1 = 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) = -AKKP1 / D 250 K=K+2 251 END IF 252 END DO 253* 254* inv(U**T) = (inv(U))**T 255* 256* inv(U**T)*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)=ONE 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**T*invD1*U11->U11 340* 341 CALL DTRMM('L','U','T','U',NNB, NNB, 342 $ ONE,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**T*invD*U01->A(CUT+I,CUT+J) 351* 352 CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, 353 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 354 355* 356* U11 = U11**T*invD1*U11 + U01**T*invD*U01 357* 358 DO I=1,NNB 359 DO J=I,NNB 360 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 361 END DO 362 END DO 363* 364* U01 = U00**T*invD0*U01 365* 366 CALL DTRMM('L',UPLO,'T','U',CUT, NNB, 367 $ ONE,A,LDA,WORK,N+NB+1) 368 369* 370* Update U01 371* 372 DO I=1,CUT 373 DO J=1,NNB 374 A(I,CUT+J)=WORK(I,J) 375 END DO 376 END DO 377* 378* Next Block 379* 380 END DO 381* 382* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 383* 384 I=1 385 DO WHILE ( I .LE. N ) 386 IF( IPIV(I) .GT. 0 ) THEN 387 IP=IPIV(I) 388 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 389 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I ) 390 ELSE 391 IP=-IPIV(I) 392 I=I+1 393 IF ( (I-1) .LT. IP) 394 $ CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) 395 IF ( (I-1) .GT. IP) 396 $ CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) 397 ENDIF 398 I=I+1 399 END DO 400 ELSE 401* 402* LOWER... 403* 404* invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 405* 406 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO ) 407* 408* inv(D) and inv(D)*inv(U) 409* 410 K=N 411 DO WHILE ( K .GE. 1 ) 412 IF( IPIV( K ).GT.0 ) THEN 413* 1 x 1 diagonal NNB 414 WORK(K,INVD) = ONE / A( K, K ) 415 WORK(K,INVD+1) = 0 416 K=K-1 417 ELSE 418* 2 x 2 diagonal NNB 419 T = WORK(K-1,1) 420 AK = A( K-1, K-1 ) / T 421 AKP1 = A( K, K ) / T 422 AKKP1 = WORK(K-1,1) / T 423 D = T*( AK*AKP1-ONE ) 424 WORK(K-1,INVD) = AKP1 / D 425 WORK(K,INVD) = AK / D 426 WORK(K,INVD+1) = -AKKP1 / D 427 WORK(K-1,INVD+1) = -AKKP1 / D 428 K=K-2 429 END IF 430 END DO 431* 432* inv(U**T) = (inv(U))**T 433* 434* inv(U**T)*inv(D)*inv(U) 435* 436 CUT=0 437 DO WHILE (CUT .LT. N) 438 NNB=NB 439 IF (CUT + NNB .GT. N) THEN 440 NNB=N-CUT 441 ELSE 442 COUNT = 0 443* count negative elements, 444 DO I=CUT+1,CUT+NNB 445 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 446 END DO 447* need a even number for a clear cut 448 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 449 END IF 450* L21 Block 451 DO I=1,N-CUT-NNB 452 DO J=1,NNB 453 WORK(I,J)=A(CUT+NNB+I,CUT+J) 454 END DO 455 END DO 456* L11 Block 457 DO I=1,NNB 458 WORK(U11+I,I)=ONE 459 DO J=I+1,NNB 460 WORK(U11+I,J)=ZERO 461 END DO 462 DO J=1,I-1 463 WORK(U11+I,J)=A(CUT+I,CUT+J) 464 END DO 465 END DO 466* 467* invD*L21 468* 469 I=N-CUT-NNB 470 DO WHILE (I .GE. 1) 471 IF (IPIV(CUT+NNB+I) > 0) THEN 472 DO J=1,NNB 473 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 474 END DO 475 I=I-1 476 ELSE 477 DO J=1,NNB 478 U01_I_J = WORK(I,J) 479 U01_IP1_J = WORK(I-1,J) 480 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 481 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 482 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 483 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 484 END DO 485 I=I-2 486 END IF 487 END DO 488* 489* invD1*L11 490* 491 I=NNB 492 DO WHILE (I .GE. 1) 493 IF (IPIV(CUT+I) > 0) THEN 494 DO J=1,NNB 495 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 496 END DO 497 I=I-1 498 ELSE 499 DO J=1,NNB 500 U11_I_J = WORK(U11+I,J) 501 U11_IP1_J = WORK(U11+I-1,J) 502 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 503 $ WORK(CUT+I,INVD+1)*U11_IP1_J 504 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 505 $ WORK(CUT+I-1,INVD)*U11_IP1_J 506 END DO 507 I=I-2 508 END IF 509 END DO 510* 511* L11**T*invD1*L11->L11 512* 513 CALL DTRMM('L',UPLO,'T','U',NNB, NNB, 514 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 515 516* 517 DO I=1,NNB 518 DO J=1,I 519 A(CUT+I,CUT+J)=WORK(U11+I,J) 520 END DO 521 END DO 522* 523 IF ( (CUT+NNB) .LT. N ) THEN 524* 525* L21**T*invD2*L21->A(CUT+I,CUT+J) 526* 527 CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) 528 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 529 530* 531* L11 = L11**T*invD1*L11 + U01**T*invD*U01 532* 533 DO I=1,NNB 534 DO J=1,I 535 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 536 END DO 537 END DO 538* 539* L01 = L22**T*invD2*L21 540* 541 CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, 542 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 543* 544* Update L21 545* 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 552 ELSE 553* 554* L11 = L11**T*invD1*L11 555* 556 DO I=1,NNB 557 DO J=1,I 558 A(CUT+I,CUT+J)=WORK(U11+I,J) 559 END DO 560 END DO 561 END IF 562* 563* Next Block 564* 565 CUT=CUT+NNB 566 END DO 567* 568* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 569* 570 I=N 571 DO WHILE ( I .GE. 1 ) 572 IF( IPIV(I) .GT. 0 ) THEN 573 IP=IPIV(I) 574 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 575 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I ) 576 ELSE 577 IP=-IPIV(I) 578 IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 579 IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I ) 580 I=I-1 581 ENDIF 582 I=I-1 583 END DO 584 END IF 585* 586 RETURN 587* 588* End of DSYTRI2X 589* 590 END 591 592