SRC/chetri2x.f

   1 *> \brief \b CHETRI2X
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CHETRI2X + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetri2x.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetri2x.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetri2x.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHETRI2X( 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            A( LDA, * ), WORK( N+NB+1,* )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CHETRI2X computes the inverse of a complex Hermitian indefinite matrix
  39 *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
  40 *> CHETRF.
  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**H;
  52 *>          = 'L':  Lower triangular, form is A = L*D*L**H.
  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 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 CHETRF.
  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 CHETRF.
  86 *> \endverbatim
  87 *>
  88 *> \param[out] WORK
  89 *> \verbatim
  90 *>          WORK is COMPLEX 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 November 2015
 117 *
 118 *> \ingroup complexHEcomputational
 119 *
 120 *  =====================================================================
 121       SUBROUTINE CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
 122 *
 123 *  -- LAPACK computational routine (version 3.6.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 2015
 127 *
 128 *     .. Scalar Arguments ..
 129       CHARACTER          UPLO
 130       INTEGER            INFO, LDA, N, NB
 131 *     ..
 132 *     .. Array Arguments ..
 133       INTEGER            IPIV( * )
 134       COMPLEX            A( LDA, * ), WORK( N+NB+1,* )
 135 *     ..
 136 *
 137 *  =====================================================================
 138 *
 139 *     .. Parameters ..
 140       REAL               ONE
 141       COMPLEX            CONE, ZERO
 142       PARAMETER          ( ONE = 1.0E+0,
 143      $                   CONE = ( 1.0E+0, 0.0E+0 ),
 144      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
 145 *     ..
 146 *     .. Local Scalars ..
 147       LOGICAL            UPPER
 148       INTEGER            I, IINFO, IP, K, CUT, NNB
 149       INTEGER            COUNT
 150       INTEGER            J, U11, INVD
 151
 152       COMPLEX   AK, AKKP1, AKP1, D, T
 153       COMPLEX   U01_I_J, U01_IP1_J
 154       COMPLEX   U11_I_J, U11_IP1_J
 155 *     ..
 156 *     .. External Functions ..
 157       LOGICAL            LSAME
 158       EXTERNAL           LSAME
 159 *     ..
 160 *     .. External Subroutines ..
 161       EXTERNAL           CSYCONV, XERBLA, CTRTRI
 162       EXTERNAL           CGEMM, CTRMM, CHESWAPR
 163 *     ..
 164 *     .. Intrinsic Functions ..
 165       INTRINSIC          MAX
 166 *     ..
 167 *     .. Executable Statements ..
 168 *
 169 *     Test the input parameters.
 170 *
 171       INFO = 0
 172       UPPER = LSAME( UPLO, 'U' )
 173       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 174          INFO = -1
 175       ELSE IF( N.LT.0 ) THEN
 176          INFO = -2
 177       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 178          INFO = -4
 179       END IF
 180 *
 181 *     Quick return if possible
 182 *
 183 *
 184       IF( INFO.NE.0 ) THEN
 185          CALL XERBLA( 'CHETRI2X', -INFO )
 186          RETURN
 187       END IF
 188       IF( N.EQ.0 )
 189      $   RETURN
 190 *
 191 *     Convert A
 192 *     Workspace got Non-diag elements of D
 193 *
 194       CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
 195 *
 196 *     Check that the diagonal matrix D is nonsingular.
 197 *
 198       IF( UPPER ) THEN
 199 *
 200 *        Upper triangular storage: examine D from bottom to top
 201 *
 202          DO INFO = N, 1, -1
 203             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 204      $         RETURN
 205          END DO
 206       ELSE
 207 *
 208 *        Lower triangular storage: examine D from top to bottom.
 209 *
 210          DO INFO = 1, N
 211             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 212      $         RETURN
 213          END DO
 214       END IF
 215       INFO = 0
 216 *
 217 *  Splitting Workspace
 218 *     U01 is a block (N,NB+1)
 219 *     The first element of U01 is in WORK(1,1)
 220 *     U11 is a block (NB+1,NB+1)
 221 *     The first element of U11 is in WORK(N+1,1)
 222       U11 = N
 223 *     INVD is a block (N,2)
 224 *     The first element of INVD is in WORK(1,INVD)
 225       INVD = NB+2
 226
 227       IF( UPPER ) THEN
 228 *
 229 *        invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
 230 *
 231         CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO )
 232 *
 233 *       inv(D) and inv(D)*inv(U)
 234 *
 235         K=1
 236         DO WHILE ( K .LE. N )
 237          IF( IPIV( K ).GT.0 ) THEN
 238 *           1 x 1 diagonal NNB
 239              WORK(K,INVD) = ONE / REAL ( A( K, K ) )
 240              WORK(K,INVD+1) = 0
 241             K=K+1
 242          ELSE
 243 *           2 x 2 diagonal NNB
 244              T = ABS ( WORK(K+1,1) )
 245              AK = REAL ( A( K, K ) ) / T
 246              AKP1 = REAL ( A( K+1, K+1 ) ) / T
 247              AKKP1 = WORK(K+1,1)  / T
 248              D = T*( AK*AKP1-ONE )
 249              WORK(K,INVD) = AKP1 / D
 250              WORK(K+1,INVD+1) = AK / D
 251              WORK(K,INVD+1) = -AKKP1 / D
 252              WORK(K+1,INVD) = CONJG (WORK(K,INVD+1) )
 253             K=K+2
 254          END IF
 255         END DO
 256 *
 257 *       inv(U**H) = (inv(U))**H
 258 *
 259 *       inv(U**H)*inv(D)*inv(U)
 260 *
 261         CUT=N
 262         DO WHILE (CUT .GT. 0)
 263            NNB=NB
 264            IF (CUT .LE. NNB) THEN
 265               NNB=CUT
 266            ELSE
 267               COUNT = 0
 268 *             count negative elements,
 269               DO I=CUT+1-NNB,CUT
 270                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
 271               END DO
 272 *             need a even number for a clear cut
 273               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
 274            END IF
 275
 276            CUT=CUT-NNB
 277 *
 278 *          U01 Block
 279 *
 280            DO I=1,CUT
 281              DO J=1,NNB
 282               WORK(I,J)=A(I,CUT+J)
 283              END DO
 284            END DO
 285 *
 286 *          U11 Block
 287 *
 288            DO I=1,NNB
 289              WORK(U11+I,I)=CONE
 290              DO J=1,I-1
 291                 WORK(U11+I,J)=ZERO
 292              END DO
 293              DO J=I+1,NNB
 294                 WORK(U11+I,J)=A(CUT+I,CUT+J)
 295              END DO
 296            END DO
 297 *
 298 *          invD*U01
 299 *
 300            I=1
 301            DO WHILE (I .LE. CUT)
 302              IF (IPIV(I) > 0) THEN
 303                 DO J=1,NNB
 304                     WORK(I,J)=WORK(I,INVD)*WORK(I,J)
 305                 END DO
 306                 I=I+1
 307              ELSE
 308                 DO J=1,NNB
 309                    U01_I_J = WORK(I,J)
 310                    U01_IP1_J = WORK(I+1,J)
 311                    WORK(I,J)=WORK(I,INVD)*U01_I_J+
 312      $                      WORK(I,INVD+1)*U01_IP1_J
 313                    WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
 314      $                      WORK(I+1,INVD+1)*U01_IP1_J
 315                 END DO
 316                 I=I+2
 317              END IF
 318            END DO
 319 *
 320 *        invD1*U11
 321 *
 322            I=1
 323            DO WHILE (I .LE. NNB)
 324              IF (IPIV(CUT+I) > 0) THEN
 325                 DO J=I,NNB
 326                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
 327                 END DO
 328                 I=I+1
 329              ELSE
 330                 DO J=I,NNB
 331                    U11_I_J = WORK(U11+I,J)
 332                    U11_IP1_J = WORK(U11+I+1,J)
 333                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
 334      $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
 335                 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
 336      $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
 337                 END DO
 338                 I=I+2
 339              END IF
 340            END DO
 341 *
 342 *       U11**H*invD1*U11->U11
 343 *
 344         CALL CTRMM('L','U','C','U',NNB, NNB,
 345      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 346 *
 347          DO I=1,NNB
 348             DO J=I,NNB
 349               A(CUT+I,CUT+J)=WORK(U11+I,J)
 350             END DO
 351          END DO
 352 *
 353 *          U01**H*invD*U01->A(CUT+I,CUT+J)
 354 *
 355          CALL CGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
 356      $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 357 *
 358 *        U11 =  U11**H*invD1*U11 + U01**H*invD*U01
 359 *
 360          DO I=1,NNB
 361             DO J=I,NNB
 362               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
 363             END DO
 364          END DO
 365 *
 366 *        U01 =  U00**H*invD0*U01
 367 *
 368          CALL CTRMM('L',UPLO,'C','U',CUT, NNB,
 369      $             CONE,A,LDA,WORK,N+NB+1)
 370
 371 *
 372 *        Update U01
 373 *
 374          DO I=1,CUT
 375            DO J=1,NNB
 376             A(I,CUT+J)=WORK(I,J)
 377            END DO
 378          END DO
 379 *
 380 *      Next Block
 381 *
 382        END DO
 383 *
 384 *        Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
 385 *
 386             I=1
 387             DO WHILE ( I .LE. N )
 388                IF( IPIV(I) .GT. 0 ) THEN
 389                   IP=IPIV(I)
 390                  IF (I .LT. IP) CALL CHESWAPR( UPLO, N, A, LDA, I ,IP )
 391                  IF (I .GT. IP) CALL CHESWAPR( UPLO, N, A, LDA, IP ,I )
 392                ELSE
 393                  IP=-IPIV(I)
 394                  I=I+1
 395                  IF ( (I-1) .LT. IP)
 396      $                  CALL CHESWAPR( UPLO, N, A, LDA, I-1 ,IP )
 397                  IF ( (I-1) .GT. IP)
 398      $                  CALL CHESWAPR( UPLO, N, A, LDA, IP ,I-1 )
 399               ENDIF
 400                I=I+1
 401             END DO
 402       ELSE
 403 *
 404 *        LOWER...
 405 *
 406 *        invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
 407 *
 408          CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO )
 409 *
 410 *       inv(D) and inv(D)*inv(U)
 411 *
 412         K=N
 413         DO WHILE ( K .GE. 1 )
 414          IF( IPIV( K ).GT.0 ) THEN
 415 *           1 x 1 diagonal NNB
 416              WORK(K,INVD) = ONE / REAL ( A( K, K ) )
 417              WORK(K,INVD+1) = 0
 418             K=K-1
 419          ELSE
 420 *           2 x 2 diagonal NNB
 421              T = ABS ( WORK(K-1,1) )
 422              AK = REAL ( A( K-1, K-1 ) ) / T
 423              AKP1 = REAL ( A( K, K ) ) / T
 424              AKKP1 = WORK(K-1,1) / T
 425              D = T*( AK*AKP1-ONE )
 426              WORK(K-1,INVD) = AKP1 / D
 427              WORK(K,INVD) = AK / D
 428              WORK(K,INVD+1) = -AKKP1 / D
 429              WORK(K-1,INVD+1) = CONJG (WORK(K,INVD+1) )
 430             K=K-2
 431          END IF
 432         END DO
 433 *
 434 *       inv(U**H) = (inv(U))**H
 435 *
 436 *       inv(U**H)*inv(D)*inv(U)
 437 *
 438         CUT=0
 439         DO WHILE (CUT .LT. N)
 440            NNB=NB
 441            IF (CUT + NNB .GE. N) THEN
 442               NNB=N-CUT
 443            ELSE
 444               COUNT = 0
 445 *             count negative elements,
 446               DO I=CUT+1,CUT+NNB
 447                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
 448               END DO
 449 *             need a even number for a clear cut
 450               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
 451            END IF
 452 *      L21 Block
 453            DO I=1,N-CUT-NNB
 454              DO J=1,NNB
 455               WORK(I,J)=A(CUT+NNB+I,CUT+J)
 456              END DO
 457            END DO
 458 *     L11 Block
 459            DO I=1,NNB
 460              WORK(U11+I,I)=CONE
 461              DO J=I+1,NNB
 462                 WORK(U11+I,J)=ZERO
 463              END DO
 464              DO J=1,I-1
 465                 WORK(U11+I,J)=A(CUT+I,CUT+J)
 466              END DO
 467            END DO
 468 *
 469 *          invD*L21
 470 *
 471            I=N-CUT-NNB
 472            DO WHILE (I .GE. 1)
 473              IF (IPIV(CUT+NNB+I) > 0) THEN
 474                 DO J=1,NNB
 475                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
 476                 END DO
 477                 I=I-1
 478              ELSE
 479                 DO J=1,NNB
 480                    U01_I_J = WORK(I,J)
 481                    U01_IP1_J = WORK(I-1,J)
 482                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
 483      $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
 484                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
 485      $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
 486                 END DO
 487                 I=I-2
 488              END IF
 489            END DO
 490 *
 491 *        invD1*L11
 492 *
 493            I=NNB
 494            DO WHILE (I .GE. 1)
 495              IF (IPIV(CUT+I) > 0) THEN
 496                 DO J=1,NNB
 497                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
 498                 END DO
 499                 I=I-1
 500              ELSE
 501                 DO J=1,NNB
 502                    U11_I_J = WORK(U11+I,J)
 503                    U11_IP1_J = WORK(U11+I-1,J)
 504                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
 505      $                      WORK(CUT+I,INVD+1)*U11_IP1_J
 506                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
 507      $                      WORK(CUT+I-1,INVD)*U11_IP1_J
 508                 END DO
 509                 I=I-2
 510              END IF
 511            END DO
 512 *
 513 *       L11**H*invD1*L11->L11
 514 *
 515         CALL CTRMM('L',UPLO,'C','U',NNB, NNB,
 516      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
 517 *
 518          DO I=1,NNB
 519             DO J=1,I
 520               A(CUT+I,CUT+J)=WORK(U11+I,J)
 521             END DO
 522          END DO
 523 *
 524         IF ( (CUT+NNB) .LT. N ) THEN
 525 *
 526 *          L21**H*invD2*L21->A(CUT+I,CUT+J)
 527 *
 528          CALL CGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
 529      $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
 530
 531 *
 532 *        L11 =  L11**H*invD1*L11 + U01**H*invD*U01
 533 *
 534          DO I=1,NNB
 535             DO J=1,I
 536               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
 537             END DO
 538          END DO
 539 *
 540 *        L01 =  L22**H*invD2*L21
 541 *
 542          CALL CTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB,
 543      $             CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
 544
 545 *      Update L21
 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        ELSE
 552 *
 553 *        L11 =  L11**H*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**H: P * inv(U**H)*inv(D)*inv(U) *P**H
 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 CHESWAPR( UPLO, N, A, LDA, I ,IP  )
 574                  IF (I .GT. IP) CALL CHESWAPR( UPLO, N, A, LDA, IP ,I )
 575                ELSE
 576                  IP=-IPIV(I)
 577                  IF ( I .LT. IP) CALL CHESWAPR( UPLO, N, A, LDA, I ,IP )
 578                  IF ( I .GT. IP) CALL CHESWAPR( 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 CHETRI2X
 588 *
 589       END
 590