SRC/zsytri2x.f

   1 *> \brief \b ZSYTRI2X
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZSYTRI2X + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri2x.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri2x.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri2x.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZSYTRI2X( 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*16         A( LDA, * ), WORK( N+NB+1,* )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> ZSYTRI2X computes the inverse of a complex symmetric indefinite matrix
  39 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
  40 *> ZSYTRF.
  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 COMPLEX*16 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 ZSYTRF.
  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 ZSYTRF.
  86 *> \endverbatim
  87 *>
  88 *> \param[out] WORK
  89 *> \verbatim
  90 *>          WORK is COMPLEX*16 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 complex16SYcomputational
 119 *
 120 *  =====================================================================
 121       SUBROUTINE ZSYTRI2X( 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       COMPLEX*16         A( LDA, * ), WORK( N+NB+1,* )
 135 *     ..
 136 *
 137 *  =====================================================================
 138 *
 139 *     .. Parameters ..
 140       COMPLEX*16         ONE, ZERO
 141       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
 142      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
 143 *     ..
 144 *     .. Local Scalars ..
 145       LOGICAL            UPPER
 146       INTEGER            I, IINFO, IP, K, CUT, NNB
 147       INTEGER            COUNT
 148       INTEGER            J, U11, INVD
 149
 150       COMPLEX*16         AK, AKKP1, AKP1, D, T
 151       COMPLEX*16         U01_I_J, U01_IP1_J
 152       COMPLEX*16         U11_I_J, U11_IP1_J
 153 *     ..
 154 *     .. External Functions ..
 155       LOGICAL            LSAME
 156       EXTERNAL           LSAME
 157 *     ..
 158 *     .. External Subroutines ..
 159       EXTERNAL           ZSYCONV, XERBLA, ZTRTRI
 160       EXTERNAL           ZGEMM, ZTRMM, ZSYSWAPR
 161 *     ..
 162 *     .. Intrinsic Functions ..
 163       INTRINSIC          MAX
 164 *     ..
 165 *     .. Executable Statements ..
 166 *
 167 *     Test the input parameters.
 168 *
 169       INFO = 0
 170       UPPER = LSAME( UPLO, 'U' )
 171       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 172          INFO = -1
 173       ELSE IF( N.LT.0 ) THEN
 174          INFO = -2
 175       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 176          INFO = -4
 177       END IF
 178 *
 179 *     Quick return if possible
 180 *
 181 *
 182       IF( INFO.NE.0 ) THEN
 183          CALL XERBLA( 'ZSYTRI2X', -INFO )
 184          RETURN
 185       END IF
 186       IF( N.EQ.0 )
 187      $   RETURN
 188 *
 189 *     Convert A
 190 *     Workspace got Non-diag elements of D
 191 *
 192       CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
 193 *
 194 *     Check that the diagonal matrix D is nonsingular.
 195 *
 196       IF( UPPER ) THEN
 197 *
 198 *        Upper triangular storage: examine D from bottom to top
 199 *
 200          DO INFO = N, 1, -1
 201             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 202      $         RETURN
 203          END DO
 204       ELSE
 205 *
 206 *        Lower triangular storage: examine D from top to bottom.
 207 *
 208          DO INFO = 1, N
 209             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 210      $         RETURN
 211          END DO
 212       END IF
 213       INFO = 0
 214 *
 215 *  Splitting Workspace
 216 *     U01 is a block (N,NB+1)
 217 *     The first element of U01 is in WORK(1,1)
 218 *     U11 is a block (NB+1,NB+1)
 219 *     The first element of U11 is in WORK(N+1,1)
 220       U11 = N
 221 *     INVD is a block (N,2)
 222 *     The first element of INVD is in WORK(1,INVD)
 223       INVD = NB+2
 224
 225       IF( UPPER ) THEN
 226 *
 227 *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
 228 *
 229         CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
 230 *
 231 *       inv(D) and inv(D)*inv(U)
 232 *
 233         K=1
 234         DO WHILE ( K .LE. N )
 235          IF( IPIV( K ).GT.0 ) THEN
 236 *           1 x 1 diagonal NNB
 237              WORK(K,INVD) = 1/  A( K, K )
 238              WORK(K,INVD+1) = 0
 239             K=K+1
 240          ELSE
 241 *           2 x 2 diagonal NNB
 242              T = WORK(K+1,1)
 243              AK = A( K, K ) / T
 244              AKP1 = A( K+1, K+1 ) / T
 245              AKKP1 = WORK(K+1,1)  / T
 246              D = T*( AK*AKP1-ONE )
 247              WORK(K,INVD) = AKP1 / D
 248              WORK(K+1,INVD+1) = AK / D
 249              WORK(K,INVD+1) = -AKKP1 / D
 250              WORK(K+1,INVD) = -AKKP1 / D
 251             K=K+2
 252          END IF
 253         END DO
 254 *
 255 *       inv(U**T) = (inv(U))**T
 256 *
 257 *       inv(U**T)*inv(D)*inv(U)
 258 *
 259         CUT=N
 260         DO WHILE (CUT .GT. 0)
 261            NNB=NB
 262            IF (CUT .LE. NNB) THEN
 263               NNB=CUT
 264            ELSE
 265               COUNT = 0
 266 *             count negative elements,
 267               DO I=CUT+1-NNB,CUT
 268                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
 269               END DO
 270 *             need a even number for a clear cut
 271               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
 272            END IF
 273
 274            CUT=CUT-NNB
 275 *
 276 *          U01 Block
 277 *
 278            DO I=1,CUT
 279              DO J=1,NNB
 280               WORK(I,J)=A(I,CUT+J)
 281              END DO
 282            END DO
 283 *
 284 *          U11 Block
 285 *
 286            DO I=1,NNB
 287              WORK(U11+I,I)=ONE
 288              DO J=1,I-1
 289                 WORK(U11+I,J)=ZERO
 290              END DO
 291              DO J=I+1,NNB
 292                 WORK(U11+I,J)=A(CUT+I,CUT+J)
 293              END DO
 294            END DO
 295 *
 296 *          invD*U01
 297 *
 298            I=1
 299            DO WHILE (I .LE. CUT)
 300              IF (IPIV(I) > 0) THEN
 301                 DO J=1,NNB
 302                     WORK(I,J)=WORK(I,INVD)*WORK(I,J)
 303                 END DO
 304                 I=I+1
 305              ELSE
 306                 DO J=1,NNB
 307                    U01_I_J = WORK(I,J)
 308                    U01_IP1_J = WORK(I+1,J)
 309                    WORK(I,J)=WORK(I,INVD)*U01_I_J+
 310      $                      WORK(I,INVD+1)*U01_IP1_J
 311                    WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
 312      $                      WORK(I+1,INVD+1)*U01_IP1_J
 313                 END DO
 314                 I=I+2
 315              END IF
 316            END DO
 317 *
 318 *        invD1*U11
 319 *
 320            I=1
 321            DO WHILE (I .LE. NNB)
 322              IF (IPIV(CUT+I) > 0) THEN
 323                 DO J=I,NNB
 324                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
 325                 END DO
 326                 I=I+1
 327              ELSE
 328                 DO J=I,NNB
 329                    U11_I_J = WORK(U11+I,J)
 330                    U11_IP1_J = WORK(U11+I+1,J)
 331                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
 332      $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
 333                 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
 334      $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
 335                 END DO
 336                 I=I+2
 337              END IF
 338            END DO
 339 *
 340 *       U11**T*invD1*U11->U11
 341 *
 342         CALL ZTRMM('L','U','T','U',NNB, NNB,
 343      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 344 *
 345          DO I=1,NNB
 346             DO J=I,NNB
 347               A(CUT+I,CUT+J)=WORK(U11+I,J)
 348             END DO
 349          END DO
 350 *
 351 *          U01**T*invD*U01->A(CUT+I,CUT+J)
 352 *
 353          CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
 354      $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 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 ZTRMM('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 ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
 389                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
 390                ELSE
 391                  IP=-IPIV(I)
 392                  I=I+1
 393                  IF ( (I-1) .LT. IP)
 394      $                  CALL ZSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
 395                  IF ( (I-1) .GT. IP)
 396      $                  CALL ZSYSWAPR( 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 ZTRTRI( 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) = 1/  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 .GE. 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 ZTRMM('L',UPLO,'T','U',NNB, NNB,
 514      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 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
 523         IF ( (CUT+NNB) .LT. N ) THEN
 524 *
 525 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
 526 *
 527          CALL ZGEMM('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 *        U01 =  L22**T*invD2*L21
 540 *
 541          CALL ZTRMM('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          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        ELSE
 551 *
 552 *        L11 =  L11**T*invD1*L11
 553 *
 554          DO I=1,NNB
 555             DO J=1,I
 556               A(CUT+I,CUT+J)=WORK(U11+I,J)
 557             END DO
 558          END DO
 559        END IF
 560 *
 561 *      Next Block
 562 *
 563            CUT=CUT+NNB
 564        END DO
 565 *
 566 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
 567 *
 568             I=N
 569             DO WHILE ( I .GE. 1 )
 570                IF( IPIV(I) .GT. 0 ) THEN
 571                   IP=IPIV(I)
 572                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP  )
 573                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
 574                ELSE
 575                  IP=-IPIV(I)
 576                  IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
 577                  IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
 578                  I=I-1
 579                ENDIF
 580                I=I-1
 581             END DO
 582       END IF
 583 *
 584       RETURN
 585 *
 586 *     End of ZSYTRI2X
 587 *
 588       END
 589