SRC/ssytri2x.f

   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+NNB+1,NNB+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 2011
 117 *
 118 *> \ingroup realSYcomputational
 119 *
 120 *  =====================================================================
 121       SUBROUTINE SSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
 122 *
 123 *  -- LAPACK computational routine (version 3.4.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 2011
 127 *
 128 *     .. Scalar Arguments ..
 129       CHARACTER          UPLO
 130       INTEGER            INFO, LDA, N, NB
 131 *     ..
 132 *     .. Array Arguments ..
 133       INTEGER            IPIV( * )
 134       REAL               A( LDA, * ), WORK( N+NB+1,* )
 135 *     ..
 136 *
 137 *  =====================================================================
 138 *
 139 *     .. Parameters ..
 140       REAL               ONE, ZERO
 141       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+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       REAL               AK, AKKP1, AKP1, D, T
 150       REAL               U01_I_J, U01_IP1_J
 151       REAL               U11_I_J, U11_IP1_J
 152 *     ..
 153 *     .. External Functions ..
 154       LOGICAL            LSAME
 155       EXTERNAL           LSAME
 156 *     ..
 157 *     .. External Subroutines ..
 158       EXTERNAL           SSYCONV, XERBLA, STRTRI
 159       EXTERNAL           SGEMM, STRMM, SSYSWAPR
 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( 'SSYTRI2X', -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 SSYCONV( 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 STRTRI( 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 STRMM('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 SGEMM('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 *        U11 =  U11**T*invD1*U11 + U01**T*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**T*invD0*U01
 364 *
 365          CALL STRMM('L',UPLO,'T','U',CUT, NNB,
 366      $             ONE,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**T: P * inv(U**T)*inv(D)*inv(U) *P**T
 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 SSYSWAPR( UPLO, N, A, LDA, I ,IP )
 388                  IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
 389                ELSE
 390                  IP=-IPIV(I)
 391                  I=I+1
 392                  IF ( (I-1) .LT. IP)
 393      $                  CALL SSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
 394                  IF ( (I-1) .GT. IP)
 395      $                  CALL SSYSWAPR( 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**T)*inv(D)*inv(U)*P**T.
 404 *
 405          CALL STRTRI( 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 /  A( K, K )
 414              WORK(K,INVD+1) = 0
 415             K=K-1
 416          ELSE
 417 *           2 x 2 diagonal NNB
 418              T = WORK(K-1,1)
 419              AK = A( K-1, K-1 ) / T
 420              AKP1 = 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) = -AKKP1 / D
 427             K=K-2
 428          END IF
 429         END DO
 430 *
 431 *       inv(U**T) = (inv(U))**T
 432 *
 433 *       inv(U**T)*inv(D)*inv(U)
 434 *
 435         CUT=0
 436         DO WHILE (CUT .LT. N)
 437            NNB=NB
 438            IF (CUT + NNB .GT. 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)=ONE
 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**T*invD1*L11->L11
 511 *
 512         CALL STRMM('L',UPLO,'T','U',NNB, NNB,
 513      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 514
 515 *
 516          DO I=1,NNB
 517             DO J=1,I
 518               A(CUT+I,CUT+J)=WORK(U11+I,J)
 519             END DO
 520          END DO
 521 *
 522         IF ( (CUT+NNB) .LT. N ) THEN
 523 *
 524 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
 525 *
 526          CALL SGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
 527      $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 528
 529 *
 530 *        L11 =  L11**T*invD1*L11 + U01**T*invD*U01
 531 *
 532          DO I=1,NNB
 533             DO J=1,I
 534               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
 535             END DO
 536          END DO
 537 *
 538 *        L01 =  L22**T*invD2*L21
 539 *
 540          CALL STRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
 541      $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
 542 *
 543 *      Update L21
 544 *
 545          DO I=1,N-CUT-NNB
 546            DO J=1,NNB
 547               A(CUT+NNB+I,CUT+J)=WORK(I,J)
 548            END DO
 549          END DO
 550
 551        ELSE
 552 *
 553 *        L11 =  L11**T*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**T: P * inv(U**T)*inv(D)*inv(U) *P**T
 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 SSYSWAPR( UPLO, N, A, LDA, I ,IP  )
 574                  IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
 575                ELSE
 576                  IP=-IPIV(I)
 577                  IF ( I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP )
 578                  IF ( I .GT. IP) CALL SSYSWAPR( 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 SSYTRI2X
 588 *
 589       END
 590