1*> \brief \b SSYTRI2X 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SSYTRI2X + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri2x.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri2x.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri2x.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SSYTRI2X( 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* REAL A( LDA, * ), WORK( N+NB+1,* ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> SSYTRI2X 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*> SSYTRF. 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 REAL 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 SSYTRF. 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 SSYTRF. 86*> \endverbatim 87*> 88*> \param[out] WORK 89*> \verbatim 90*> WORK is REAL 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 realSYcomputational 117* 118* ===================================================================== 119 SUBROUTINE SSYTRI2X( 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 REAL A( LDA, * ), WORK( N+NB+1,* ) 132* .. 133* 134* ===================================================================== 135* 136* .. Parameters .. 137 REAL ONE, ZERO 138 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+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 REAL AK, AKKP1, AKP1, D, T 147 REAL U01_I_J, U01_IP1_J 148 REAL U11_I_J, U11_IP1_J 149* .. 150* .. External Functions .. 151 LOGICAL LSAME 152 EXTERNAL LSAME 153* .. 154* .. External Subroutines .. 155 EXTERNAL SSYCONV, XERBLA, STRTRI 156 EXTERNAL SGEMM, STRMM, SSYSWAPR 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( 'SSYTRI2X', -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 SSYCONV( 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 STRTRI( 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 STRMM('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 SGEMM('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* U11 = U11**T*invD1*U11 + U01**T*invD*U01 353* 354 DO I=1,NNB 355 DO J=I,NNB 356 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 357 END DO 358 END DO 359* 360* U01 = U00**T*invD0*U01 361* 362 CALL STRMM('L',UPLO,'T','U',CUT, NNB, 363 $ ONE,A,LDA,WORK,N+NB+1) 364 365* 366* Update U01 367* 368 DO I=1,CUT 369 DO J=1,NNB 370 A(I,CUT+J)=WORK(I,J) 371 END DO 372 END DO 373* 374* Next Block 375* 376 END DO 377* 378* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 379* 380 I=1 381 DO WHILE ( I .LE. N ) 382 IF( IPIV(I) .GT. 0 ) THEN 383 IP=IPIV(I) 384 IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP ) 385 IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I ) 386 ELSE 387 IP=-IPIV(I) 388 I=I+1 389 IF ( (I-1) .LT. IP) 390 $ CALL SSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) 391 IF ( (I-1) .GT. IP) 392 $ CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) 393 ENDIF 394 I=I+1 395 END DO 396 ELSE 397* 398* LOWER... 399* 400* invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 401* 402 CALL STRTRI( UPLO, 'U', N, A, LDA, INFO ) 403* 404* inv(D) and inv(D)*inv(U) 405* 406 K=N 407 DO WHILE ( K .GE. 1 ) 408 IF( IPIV( K ).GT.0 ) THEN 409* 1 x 1 diagonal NNB 410 WORK(K,INVD) = ONE / A( K, K ) 411 WORK(K,INVD+1) = 0 412 K=K-1 413 ELSE 414* 2 x 2 diagonal NNB 415 T = WORK(K-1,1) 416 AK = A( K-1, K-1 ) / T 417 AKP1 = A( K, K ) / T 418 AKKP1 = WORK(K-1,1) / T 419 D = T*( AK*AKP1-ONE ) 420 WORK(K-1,INVD) = AKP1 / D 421 WORK(K,INVD) = AK / D 422 WORK(K,INVD+1) = -AKKP1 / D 423 WORK(K-1,INVD+1) = -AKKP1 / D 424 K=K-2 425 END IF 426 END DO 427* 428* inv(U**T) = (inv(U))**T 429* 430* inv(U**T)*inv(D)*inv(U) 431* 432 CUT=0 433 DO WHILE (CUT .LT. N) 434 NNB=NB 435 IF (CUT + NNB .GT. N) THEN 436 NNB=N-CUT 437 ELSE 438 COUNT = 0 439* count negative elements, 440 DO I=CUT+1,CUT+NNB 441 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 442 END DO 443* need a even number for a clear cut 444 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 445 END IF 446* L21 Block 447 DO I=1,N-CUT-NNB 448 DO J=1,NNB 449 WORK(I,J)=A(CUT+NNB+I,CUT+J) 450 END DO 451 END DO 452* L11 Block 453 DO I=1,NNB 454 WORK(U11+I,I)=ONE 455 DO J=I+1,NNB 456 WORK(U11+I,J)=ZERO 457 END DO 458 DO J=1,I-1 459 WORK(U11+I,J)=A(CUT+I,CUT+J) 460 END DO 461 END DO 462* 463* invD*L21 464* 465 I=N-CUT-NNB 466 DO WHILE (I .GE. 1) 467 IF (IPIV(CUT+NNB+I) > 0) THEN 468 DO J=1,NNB 469 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 470 END DO 471 I=I-1 472 ELSE 473 DO J=1,NNB 474 U01_I_J = WORK(I,J) 475 U01_IP1_J = WORK(I-1,J) 476 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 477 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 478 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 479 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 480 END DO 481 I=I-2 482 END IF 483 END DO 484* 485* invD1*L11 486* 487 I=NNB 488 DO WHILE (I .GE. 1) 489 IF (IPIV(CUT+I) > 0) THEN 490 DO J=1,NNB 491 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 492 END DO 493 I=I-1 494 ELSE 495 DO J=1,NNB 496 U11_I_J = WORK(U11+I,J) 497 U11_IP1_J = WORK(U11+I-1,J) 498 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 499 $ WORK(CUT+I,INVD+1)*U11_IP1_J 500 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 501 $ WORK(CUT+I-1,INVD)*U11_IP1_J 502 END DO 503 I=I-2 504 END IF 505 END DO 506* 507* L11**T*invD1*L11->L11 508* 509 CALL STRMM('L',UPLO,'T','U',NNB, NNB, 510 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 511 512* 513 DO I=1,NNB 514 DO J=1,I 515 A(CUT+I,CUT+J)=WORK(U11+I,J) 516 END DO 517 END DO 518* 519 IF ( (CUT+NNB) .LT. N ) THEN 520* 521* L21**T*invD2*L21->A(CUT+I,CUT+J) 522* 523 CALL SGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) 524 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 525 526* 527* L11 = L11**T*invD1*L11 + U01**T*invD*U01 528* 529 DO I=1,NNB 530 DO J=1,I 531 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 532 END DO 533 END DO 534* 535* L01 = L22**T*invD2*L21 536* 537 CALL STRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, 538 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 539* 540* Update L21 541* 542 DO I=1,N-CUT-NNB 543 DO J=1,NNB 544 A(CUT+NNB+I,CUT+J)=WORK(I,J) 545 END DO 546 END DO 547 548 ELSE 549* 550* L11 = L11**T*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**T: P * inv(U**T)*inv(D)*inv(U) *P**T 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 SSYSWAPR( UPLO, N, A, LDA, I ,IP ) 571 IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I ) 572 ELSE 573 IP=-IPIV(I) 574 IF ( I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP ) 575 IF ( I .GT. IP) CALL SSYSWAPR( 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 SSYTRI2X 585* 586 END 587 588