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*> \ingroup doubleSYcomputational 117* 118* ===================================================================== 119 SUBROUTINE DSYTRI2X( 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 DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* ) 132* .. 133* 134* ===================================================================== 135* 136* .. Parameters .. 137 DOUBLE PRECISION ONE, ZERO 138 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 139* .. 140* .. Local Scalars .. 141 LOGICAL UPPER 142 INTEGER I, IINFO, IP, K, CUT, NNB 143 INTEGER COUNT 144 INTEGER J, U11, INVD 145 146 DOUBLE PRECISION AK, AKKP1, AKP1, D, T 147 DOUBLE PRECISION U01_I_J, U01_IP1_J 148 DOUBLE PRECISION U11_I_J, U11_IP1_J 149* .. 150* .. External Functions .. 151 LOGICAL LSAME 152 EXTERNAL LSAME 153* .. 154* .. External Subroutines .. 155 EXTERNAL DSYCONV, XERBLA, DTRTRI 156 EXTERNAL DGEMM, DTRMM, DSYSWAPR 157* .. 158* .. Intrinsic Functions .. 159 INTRINSIC MAX 160* .. 161* .. Executable Statements .. 162* 163* Test the input parameters. 164* 165 INFO = 0 166 UPPER = LSAME( UPLO, 'U' ) 167 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 168 INFO = -1 169 ELSE IF( N.LT.0 ) THEN 170 INFO = -2 171 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 172 INFO = -4 173 END IF 174* 175* Quick return if possible 176* 177* 178 IF( INFO.NE.0 ) THEN 179 CALL XERBLA( 'DSYTRI2X', -INFO ) 180 RETURN 181 END IF 182 IF( N.EQ.0 ) 183 $ RETURN 184* 185* Convert A 186* Workspace got Non-diag elements of D 187* 188 CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 189* 190* Check that the diagonal matrix D is nonsingular. 191* 192 IF( UPPER ) THEN 193* 194* Upper triangular storage: examine D from bottom to top 195* 196 DO INFO = N, 1, -1 197 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 198 $ RETURN 199 END DO 200 ELSE 201* 202* Lower triangular storage: examine D from top to bottom. 203* 204 DO INFO = 1, N 205 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 206 $ RETURN 207 END DO 208 END IF 209 INFO = 0 210* 211* Splitting Workspace 212* U01 is a block (N,NB+1) 213* The first element of U01 is in WORK(1,1) 214* U11 is a block (NB+1,NB+1) 215* The first element of U11 is in WORK(N+1,1) 216 U11 = N 217* INVD is a block (N,2) 218* The first element of INVD is in WORK(1,INVD) 219 INVD = NB+2 220 221 IF( UPPER ) THEN 222* 223* invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 224* 225 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO ) 226* 227* inv(D) and inv(D)*inv(U) 228* 229 K=1 230 DO WHILE ( K .LE. N ) 231 IF( IPIV( K ).GT.0 ) THEN 232* 1 x 1 diagonal NNB 233 WORK(K,INVD) = ONE / A( K, K ) 234 WORK(K,INVD+1) = 0 235 K=K+1 236 ELSE 237* 2 x 2 diagonal NNB 238 T = WORK(K+1,1) 239 AK = A( K, K ) / T 240 AKP1 = A( K+1, K+1 ) / T 241 AKKP1 = WORK(K+1,1) / T 242 D = T*( AK*AKP1-ONE ) 243 WORK(K,INVD) = AKP1 / D 244 WORK(K+1,INVD+1) = AK / D 245 WORK(K,INVD+1) = -AKKP1 / D 246 WORK(K+1,INVD) = -AKKP1 / D 247 K=K+2 248 END IF 249 END DO 250* 251* inv(U**T) = (inv(U))**T 252* 253* inv(U**T)*inv(D)*inv(U) 254* 255 CUT=N 256 DO WHILE (CUT .GT. 0) 257 NNB=NB 258 IF (CUT .LE. NNB) THEN 259 NNB=CUT 260 ELSE 261 COUNT = 0 262* count negative elements, 263 DO I=CUT+1-NNB,CUT 264 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 265 END DO 266* need a even number for a clear cut 267 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 268 END IF 269 270 CUT=CUT-NNB 271* 272* U01 Block 273* 274 DO I=1,CUT 275 DO J=1,NNB 276 WORK(I,J)=A(I,CUT+J) 277 END DO 278 END DO 279* 280* U11 Block 281* 282 DO I=1,NNB 283 WORK(U11+I,I)=ONE 284 DO J=1,I-1 285 WORK(U11+I,J)=ZERO 286 END DO 287 DO J=I+1,NNB 288 WORK(U11+I,J)=A(CUT+I,CUT+J) 289 END DO 290 END DO 291* 292* invD*U01 293* 294 I=1 295 DO WHILE (I .LE. CUT) 296 IF (IPIV(I) > 0) THEN 297 DO J=1,NNB 298 WORK(I,J)=WORK(I,INVD)*WORK(I,J) 299 END DO 300 I=I+1 301 ELSE 302 DO J=1,NNB 303 U01_I_J = WORK(I,J) 304 U01_IP1_J = WORK(I+1,J) 305 WORK(I,J)=WORK(I,INVD)*U01_I_J+ 306 $ WORK(I,INVD+1)*U01_IP1_J 307 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+ 308 $ WORK(I+1,INVD+1)*U01_IP1_J 309 END DO 310 I=I+2 311 END IF 312 END DO 313* 314* invD1*U11 315* 316 I=1 317 DO WHILE (I .LE. NNB) 318 IF (IPIV(CUT+I) > 0) THEN 319 DO J=I,NNB 320 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 321 END DO 322 I=I+1 323 ELSE 324 DO J=I,NNB 325 U11_I_J = WORK(U11+I,J) 326 U11_IP1_J = WORK(U11+I+1,J) 327 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 328 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J) 329 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+ 330 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J 331 END DO 332 I=I+2 333 END IF 334 END DO 335* 336* U11**T*invD1*U11->U11 337* 338 CALL DTRMM('L','U','T','U',NNB, NNB, 339 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 340* 341 DO I=1,NNB 342 DO J=I,NNB 343 A(CUT+I,CUT+J)=WORK(U11+I,J) 344 END DO 345 END DO 346* 347* U01**T*invD*U01->A(CUT+I,CUT+J) 348* 349 CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, 350 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 351 352* 353* U11 = U11**T*invD1*U11 + U01**T*invD*U01 354* 355 DO I=1,NNB 356 DO J=I,NNB 357 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 358 END DO 359 END DO 360* 361* U01 = U00**T*invD0*U01 362* 363 CALL DTRMM('L',UPLO,'T','U',CUT, NNB, 364 $ ONE,A,LDA,WORK,N+NB+1) 365 366* 367* Update U01 368* 369 DO I=1,CUT 370 DO J=1,NNB 371 A(I,CUT+J)=WORK(I,J) 372 END DO 373 END DO 374* 375* Next Block 376* 377 END DO 378* 379* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 380* 381 I=1 382 DO WHILE ( I .LE. N ) 383 IF( IPIV(I) .GT. 0 ) THEN 384 IP=IPIV(I) 385 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 386 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I ) 387 ELSE 388 IP=-IPIV(I) 389 I=I+1 390 IF ( (I-1) .LT. IP) 391 $ CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) 392 IF ( (I-1) .GT. IP) 393 $ CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) 394 ENDIF 395 I=I+1 396 END DO 397 ELSE 398* 399* LOWER... 400* 401* invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 402* 403 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO ) 404* 405* inv(D) and inv(D)*inv(U) 406* 407 K=N 408 DO WHILE ( K .GE. 1 ) 409 IF( IPIV( K ).GT.0 ) THEN 410* 1 x 1 diagonal NNB 411 WORK(K,INVD) = ONE / A( K, K ) 412 WORK(K,INVD+1) = 0 413 K=K-1 414 ELSE 415* 2 x 2 diagonal NNB 416 T = WORK(K-1,1) 417 AK = A( K-1, K-1 ) / T 418 AKP1 = A( K, K ) / T 419 AKKP1 = WORK(K-1,1) / T 420 D = T*( AK*AKP1-ONE ) 421 WORK(K-1,INVD) = AKP1 / D 422 WORK(K,INVD) = AK / D 423 WORK(K,INVD+1) = -AKKP1 / D 424 WORK(K-1,INVD+1) = -AKKP1 / D 425 K=K-2 426 END IF 427 END DO 428* 429* inv(U**T) = (inv(U))**T 430* 431* inv(U**T)*inv(D)*inv(U) 432* 433 CUT=0 434 DO WHILE (CUT .LT. N) 435 NNB=NB 436 IF (CUT + NNB .GT. N) THEN 437 NNB=N-CUT 438 ELSE 439 COUNT = 0 440* count negative elements, 441 DO I=CUT+1,CUT+NNB 442 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 443 END DO 444* need a even number for a clear cut 445 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 446 END IF 447* L21 Block 448 DO I=1,N-CUT-NNB 449 DO J=1,NNB 450 WORK(I,J)=A(CUT+NNB+I,CUT+J) 451 END DO 452 END DO 453* L11 Block 454 DO I=1,NNB 455 WORK(U11+I,I)=ONE 456 DO J=I+1,NNB 457 WORK(U11+I,J)=ZERO 458 END DO 459 DO J=1,I-1 460 WORK(U11+I,J)=A(CUT+I,CUT+J) 461 END DO 462 END DO 463* 464* invD*L21 465* 466 I=N-CUT-NNB 467 DO WHILE (I .GE. 1) 468 IF (IPIV(CUT+NNB+I) > 0) THEN 469 DO J=1,NNB 470 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 471 END DO 472 I=I-1 473 ELSE 474 DO J=1,NNB 475 U01_I_J = WORK(I,J) 476 U01_IP1_J = WORK(I-1,J) 477 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 478 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 479 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 480 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 481 END DO 482 I=I-2 483 END IF 484 END DO 485* 486* invD1*L11 487* 488 I=NNB 489 DO WHILE (I .GE. 1) 490 IF (IPIV(CUT+I) > 0) THEN 491 DO J=1,NNB 492 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 493 END DO 494 I=I-1 495 ELSE 496 DO J=1,NNB 497 U11_I_J = WORK(U11+I,J) 498 U11_IP1_J = WORK(U11+I-1,J) 499 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 500 $ WORK(CUT+I,INVD+1)*U11_IP1_J 501 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 502 $ WORK(CUT+I-1,INVD)*U11_IP1_J 503 END DO 504 I=I-2 505 END IF 506 END DO 507* 508* L11**T*invD1*L11->L11 509* 510 CALL DTRMM('L',UPLO,'T','U',NNB, NNB, 511 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 512 513* 514 DO I=1,NNB 515 DO J=1,I 516 A(CUT+I,CUT+J)=WORK(U11+I,J) 517 END DO 518 END DO 519* 520 IF ( (CUT+NNB) .LT. N ) THEN 521* 522* L21**T*invD2*L21->A(CUT+I,CUT+J) 523* 524 CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) 525 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 526 527* 528* L11 = L11**T*invD1*L11 + U01**T*invD*U01 529* 530 DO I=1,NNB 531 DO J=1,I 532 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 533 END DO 534 END DO 535* 536* L01 = L22**T*invD2*L21 537* 538 CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, 539 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 540* 541* Update L21 542* 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 549 ELSE 550* 551* L11 = L11**T*invD1*L11 552* 553 DO I=1,NNB 554 DO J=1,I 555 A(CUT+I,CUT+J)=WORK(U11+I,J) 556 END DO 557 END DO 558 END IF 559* 560* Next Block 561* 562 CUT=CUT+NNB 563 END DO 564* 565* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 566* 567 I=N 568 DO WHILE ( I .GE. 1 ) 569 IF( IPIV(I) .GT. 0 ) THEN 570 IP=IPIV(I) 571 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 572 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I ) 573 ELSE 574 IP=-IPIV(I) 575 IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 576 IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I ) 577 I=I-1 578 ENDIF 579 I=I-1 580 END DO 581 END IF 582* 583 RETURN 584* 585* End of DSYTRI2X 586* 587 END 588 589