LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
chetri2x.f
Go to the documentation of this file.
00001 *> \brief \b CHETRI2X
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CHETRI2X + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetri2x.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetri2x.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetri2x.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, LDA, N, NB
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IPIV( * )
00029 *       COMPLEX            A( LDA, * ), WORK( N+NB+1,* )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> CHETRI2X computes the inverse of a complex Hermitian indefinite matrix
00039 *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
00040 *> CHETRF.
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] UPLO
00047 *> \verbatim
00048 *>          UPLO is CHARACTER*1
00049 *>          Specifies whether the details of the factorization are stored
00050 *>          as an upper or lower triangular matrix.
00051 *>          = 'U':  Upper triangular, form is A = U*D*U**H;
00052 *>          = 'L':  Lower triangular, form is A = L*D*L**H.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] N
00056 *> \verbatim
00057 *>          N is INTEGER
00058 *>          The order of the matrix A.  N >= 0.
00059 *> \endverbatim
00060 *>
00061 *> \param[in,out] A
00062 *> \verbatim
00063 *>          A is COMPLEX array, dimension (LDA,N)
00064 *>          On entry, the NNB diagonal matrix D and the multipliers
00065 *>          used to obtain the factor U or L as computed by CHETRF.
00066 *>
00067 *>          On exit, if INFO = 0, the (symmetric) inverse of the original
00068 *>          matrix.  If UPLO = 'U', the upper triangular part of the
00069 *>          inverse is formed and the part of A below the diagonal is not
00070 *>          referenced; if UPLO = 'L' the lower triangular part of the
00071 *>          inverse is formed and the part of A above the diagonal is
00072 *>          not referenced.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>          The leading dimension of the array A.  LDA >= max(1,N).
00079 *> \endverbatim
00080 *>
00081 *> \param[in] IPIV
00082 *> \verbatim
00083 *>          IPIV is INTEGER array, dimension (N)
00084 *>          Details of the interchanges and the NNB structure of D
00085 *>          as determined by CHETRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[out] WORK
00089 *> \verbatim
00090 *>          WORK is COMPLEX array, dimension (N+NNB+1,NNB+3)
00091 *> \endverbatim
00092 *>
00093 *> \param[in] NB
00094 *> \verbatim
00095 *>          NB is INTEGER
00096 *>          Block size
00097 *> \endverbatim
00098 *>
00099 *> \param[out] INFO
00100 *> \verbatim
00101 *>          INFO is INTEGER
00102 *>          = 0: successful exit
00103 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00104 *>          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
00105 *>               inverse could not be computed.
00106 *> \endverbatim
00107 *
00108 *  Authors:
00109 *  ========
00110 *
00111 *> \author Univ. of Tennessee 
00112 *> \author Univ. of California Berkeley 
00113 *> \author Univ. of Colorado Denver 
00114 *> \author NAG Ltd. 
00115 *
00116 *> \date November 2011
00117 *
00118 *> \ingroup complexHEcomputational
00119 *
00120 *  =====================================================================
00121       SUBROUTINE CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
00122 *
00123 *  -- LAPACK computational routine (version 3.4.0) --
00124 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00125 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00126 *     November 2011
00127 *
00128 *     .. Scalar Arguments ..
00129       CHARACTER          UPLO
00130       INTEGER            INFO, LDA, N, NB
00131 *     ..
00132 *     .. Array Arguments ..
00133       INTEGER            IPIV( * )
00134       COMPLEX            A( LDA, * ), WORK( N+NB+1,* )
00135 *     ..
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       REAL               ONE
00141       COMPLEX            CONE, ZERO
00142       PARAMETER          ( ONE = 1.0E+0,
00143      $                   CONE = ( 1.0E+0, 0.0E+0 ),
00144      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
00145 *     ..
00146 *     .. Local Scalars ..
00147       LOGICAL            UPPER
00148       INTEGER            I, IINFO, IP, K, CUT, NNB
00149       INTEGER            COUNT
00150       INTEGER            J, U11, INVD
00151 
00152       COMPLEX   AK, AKKP1, AKP1, D, T
00153       COMPLEX   U01_I_J, U01_IP1_J
00154       COMPLEX   U11_I_J, U11_IP1_J
00155 *     ..
00156 *     .. External Functions ..
00157       LOGICAL            LSAME
00158       EXTERNAL           LSAME
00159 *     ..
00160 *     .. External Subroutines ..
00161       EXTERNAL           CSYCONV, XERBLA, CTRTRI
00162       EXTERNAL           CGEMM, CTRMM, CHESWAPR
00163 *     ..
00164 *     .. Intrinsic Functions ..
00165       INTRINSIC          MAX
00166 *     ..
00167 *     .. Executable Statements ..
00168 *
00169 *     Test the input parameters.
00170 *
00171       INFO = 0
00172       UPPER = LSAME( UPLO, 'U' )
00173       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00174          INFO = -1
00175       ELSE IF( N.LT.0 ) THEN
00176          INFO = -2
00177       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00178          INFO = -4
00179       END IF
00180 *
00181 *     Quick return if possible
00182 *
00183 *
00184       IF( INFO.NE.0 ) THEN
00185          CALL XERBLA( 'CHETRI2X', -INFO )
00186          RETURN
00187       END IF
00188       IF( N.EQ.0 )
00189      $   RETURN
00190 *
00191 *     Convert A
00192 *     Workspace got Non-diag elements of D
00193 *
00194       CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00195 *
00196 *     Check that the diagonal matrix D is nonsingular.
00197 *
00198       IF( UPPER ) THEN
00199 *
00200 *        Upper triangular storage: examine D from bottom to top
00201 *
00202          DO INFO = N, 1, -1
00203             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00204      $         RETURN
00205          END DO
00206       ELSE
00207 *
00208 *        Lower triangular storage: examine D from top to bottom.
00209 *
00210          DO INFO = 1, N
00211             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00212      $         RETURN
00213          END DO
00214       END IF
00215       INFO = 0
00216 *
00217 *  Splitting Workspace
00218 *     U01 is a block (N,NB+1) 
00219 *     The first element of U01 is in WORK(1,1)
00220 *     U11 is a block (NB+1,NB+1)
00221 *     The first element of U11 is in WORK(N+1,1)
00222       U11 = N
00223 *     INVD is a block (N,2)
00224 *     The first element of INVD is in WORK(1,INVD)
00225       INVD = NB+2
00226 
00227       IF( UPPER ) THEN
00228 *
00229 *        invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
00230 *
00231         CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO )
00232 *
00233 *       inv(D) and inv(D)*inv(U)
00234 * 
00235         K=1
00236         DO WHILE ( K .LE. N )
00237          IF( IPIV( K ).GT.0 ) THEN
00238 *           1 x 1 diagonal NNB
00239              WORK(K,INVD) = ONE / REAL ( A( K, K ) )
00240              WORK(K,INVD+1) = 0
00241             K=K+1
00242          ELSE
00243 *           2 x 2 diagonal NNB
00244              T = ABS ( WORK(K+1,1) )
00245              AK = REAL ( A( K, K ) ) / T
00246              AKP1 = REAL ( A( K+1, K+1 ) ) / T
00247              AKKP1 = WORK(K+1,1)  / T
00248              D = T*( AK*AKP1-ONE )
00249              WORK(K,INVD) = AKP1 / D
00250              WORK(K+1,INVD+1) = AK / D
00251              WORK(K,INVD+1) = -AKKP1 / D  
00252              WORK(K+1,INVD) = CONJG (WORK(K,INVD+1) )
00253             K=K+2
00254          END IF
00255         END DO
00256 *
00257 *       inv(U**H) = (inv(U))**H
00258 *
00259 *       inv(U**H)*inv(D)*inv(U)
00260 *
00261         CUT=N
00262         DO WHILE (CUT .GT. 0)
00263            NNB=NB
00264            IF (CUT .LE. NNB) THEN
00265               NNB=CUT
00266            ELSE
00267               COUNT = 0
00268 *             count negative elements, 
00269               DO I=CUT+1-NNB,CUT
00270                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00271               END DO
00272 *             need a even number for a clear cut
00273               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00274            END IF
00275 
00276            CUT=CUT-NNB
00277 *
00278 *          U01 Block 
00279 *
00280            DO I=1,CUT
00281              DO J=1,NNB
00282               WORK(I,J)=A(I,CUT+J)
00283              END DO
00284            END DO
00285 *
00286 *          U11 Block
00287 *
00288            DO I=1,NNB
00289              WORK(U11+I,I)=CONE
00290              DO J=1,I-1
00291                 WORK(U11+I,J)=ZERO
00292              END DO
00293              DO J=I+1,NNB
00294                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00295              END DO
00296            END DO
00297 *
00298 *          invD*U01
00299 *
00300            I=1
00301            DO WHILE (I .LE. CUT)
00302              IF (IPIV(I) > 0) THEN
00303                 DO J=1,NNB
00304                     WORK(I,J)=WORK(I,INVD)*WORK(I,J)
00305                 END DO
00306                 I=I+1
00307              ELSE
00308                 DO J=1,NNB
00309                    U01_I_J = WORK(I,J)
00310                    U01_IP1_J = WORK(I+1,J)
00311                    WORK(I,J)=WORK(I,INVD)*U01_I_J+
00312      $                      WORK(I,INVD+1)*U01_IP1_J
00313                    WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
00314      $                      WORK(I+1,INVD+1)*U01_IP1_J
00315                 END DO
00316                 I=I+2
00317              END IF
00318            END DO
00319 *
00320 *        invD1*U11
00321 *
00322            I=1
00323            DO WHILE (I .LE. NNB)
00324              IF (IPIV(CUT+I) > 0) THEN
00325                 DO J=I,NNB
00326                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00327                 END DO
00328                 I=I+1
00329              ELSE
00330                 DO J=I,NNB
00331                    U11_I_J = WORK(U11+I,J)
00332                    U11_IP1_J = WORK(U11+I+1,J)
00333                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00334      $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
00335                 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
00336      $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
00337                 END DO
00338                 I=I+2
00339              END IF
00340            END DO
00341 *    
00342 *       U11**H*invD1*U11->U11
00343 *
00344         CALL CTRMM('L','U','C','U',NNB, NNB,
00345      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00346 *
00347          DO I=1,NNB
00348             DO J=I,NNB
00349               A(CUT+I,CUT+J)=WORK(U11+I,J)
00350             END DO
00351          END DO
00352 *
00353 *          U01**H*invD*U01->A(CUT+I,CUT+J)
00354 *
00355          CALL CGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
00356      $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00357 *
00358 *        U11 =  U11**H*invD1*U11 + U01**H*invD*U01
00359 *
00360          DO I=1,NNB
00361             DO J=I,NNB
00362               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00363             END DO
00364          END DO
00365 *
00366 *        U01 =  U00**H*invD0*U01
00367 *
00368          CALL CTRMM('L',UPLO,'C','U',CUT, NNB,
00369      $             CONE,A,LDA,WORK,N+NB+1)
00370 
00371 *
00372 *        Update U01
00373 *
00374          DO I=1,CUT
00375            DO J=1,NNB
00376             A(I,CUT+J)=WORK(I,J)
00377            END DO
00378          END DO
00379 *
00380 *      Next Block
00381 *
00382        END DO
00383 *
00384 *        Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
00385 *  
00386             I=1
00387             DO WHILE ( I .LE. N )
00388                IF( IPIV(I) .GT. 0 ) THEN
00389                   IP=IPIV(I)
00390                  IF (I .LT. IP) CALL CHESWAPR( UPLO, N, A, LDA, I ,IP )
00391                  IF (I .GT. IP) CALL CHESWAPR( UPLO, N, A, LDA, IP ,I )
00392                ELSE
00393                  IP=-IPIV(I)
00394                  I=I+1
00395                  IF ( (I-1) .LT. IP) 
00396      $                  CALL CHESWAPR( UPLO, N, A, LDA, I-1 ,IP )
00397                  IF ( (I-1) .GT. IP) 
00398      $                  CALL CHESWAPR( UPLO, N, A, LDA, IP ,I-1 )
00399               ENDIF
00400                I=I+1
00401             END DO
00402       ELSE
00403 *
00404 *        LOWER...
00405 *
00406 *        invA = P * inv(U**H)*inv(D)*inv(U)*P**H.
00407 *
00408          CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO )
00409 *
00410 *       inv(D) and inv(D)*inv(U)
00411 * 
00412         K=N
00413         DO WHILE ( K .GE. 1 )
00414          IF( IPIV( K ).GT.0 ) THEN
00415 *           1 x 1 diagonal NNB
00416              WORK(K,INVD) = ONE / REAL ( A( K, K ) )
00417              WORK(K,INVD+1) = 0
00418             K=K-1
00419          ELSE
00420 *           2 x 2 diagonal NNB
00421              T = ABS ( WORK(K-1,1) )
00422              AK = REAL ( A( K-1, K-1 ) ) / T
00423              AKP1 = REAL ( A( K, K ) ) / T
00424              AKKP1 = WORK(K-1,1) / T
00425              D = T*( AK*AKP1-ONE )
00426              WORK(K-1,INVD) = AKP1 / D
00427              WORK(K,INVD) = AK / D
00428              WORK(K,INVD+1) = -AKKP1 / D  
00429              WORK(K-1,INVD+1) = CONJG (WORK(K,INVD+1) )
00430             K=K-2
00431          END IF
00432         END DO
00433 *
00434 *       inv(U**H) = (inv(U))**H
00435 *
00436 *       inv(U**H)*inv(D)*inv(U)
00437 *
00438         CUT=0
00439         DO WHILE (CUT .LT. N)
00440            NNB=NB
00441            IF (CUT + NNB .GE. N) THEN
00442               NNB=N-CUT
00443            ELSE
00444               COUNT = 0
00445 *             count negative elements, 
00446               DO I=CUT+1,CUT+NNB
00447                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00448               END DO
00449 *             need a even number for a clear cut
00450               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00451            END IF
00452 *      L21 Block
00453            DO I=1,N-CUT-NNB
00454              DO J=1,NNB
00455               WORK(I,J)=A(CUT+NNB+I,CUT+J)
00456              END DO
00457            END DO
00458 *     L11 Block
00459            DO I=1,NNB
00460              WORK(U11+I,I)=CONE
00461              DO J=I+1,NNB
00462                 WORK(U11+I,J)=ZERO
00463              END DO
00464              DO J=1,I-1
00465                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00466              END DO
00467            END DO
00468 *
00469 *          invD*L21
00470 *
00471            I=N-CUT-NNB
00472            DO WHILE (I .GE. 1)
00473              IF (IPIV(CUT+NNB+I) > 0) THEN
00474                 DO J=1,NNB
00475                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
00476                 END DO
00477                 I=I-1
00478              ELSE
00479                 DO J=1,NNB
00480                    U01_I_J = WORK(I,J)
00481                    U01_IP1_J = WORK(I-1,J)
00482                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
00483      $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
00484                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
00485      $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
00486                 END DO
00487                 I=I-2
00488              END IF
00489            END DO
00490 *
00491 *        invD1*L11
00492 *
00493            I=NNB
00494            DO WHILE (I .GE. 1)
00495              IF (IPIV(CUT+I) > 0) THEN
00496                 DO J=1,NNB
00497                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00498                 END DO
00499                 I=I-1
00500              ELSE
00501                 DO J=1,NNB
00502                    U11_I_J = WORK(U11+I,J)
00503                    U11_IP1_J = WORK(U11+I-1,J)
00504                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00505      $                      WORK(CUT+I,INVD+1)*U11_IP1_J
00506                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
00507      $                      WORK(CUT+I-1,INVD)*U11_IP1_J
00508                 END DO
00509                 I=I-2
00510              END IF
00511            END DO
00512 *    
00513 *       L11**H*invD1*L11->L11
00514 *
00515         CALL CTRMM('L',UPLO,'C','U',NNB, NNB,
00516      $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00517 *
00518          DO I=1,NNB
00519             DO J=1,I
00520               A(CUT+I,CUT+J)=WORK(U11+I,J)
00521             END DO
00522          END DO
00523 *
00524         IF ( (CUT+NNB) .LT. N ) THEN
00525 *
00526 *          L21**H*invD2*L21->A(CUT+I,CUT+J)
00527 *
00528          CALL CGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
00529      $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00530        
00531 *
00532 *        L11 =  L11**H*invD1*L11 + U01**H*invD*U01
00533 *
00534          DO I=1,NNB
00535             DO J=1,I
00536               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00537             END DO
00538          END DO
00539 *
00540 *        L01 =  L22**H*invD2*L21
00541 *
00542          CALL CTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB,
00543      $             CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
00544 
00545 *      Update L21
00546          DO I=1,N-CUT-NNB
00547            DO J=1,NNB
00548               A(CUT+NNB+I,CUT+J)=WORK(I,J)
00549            END DO
00550          END DO
00551        ELSE
00552 *
00553 *        L11 =  L11**H*invD1*L11
00554 *
00555          DO I=1,NNB
00556             DO J=1,I
00557               A(CUT+I,CUT+J)=WORK(U11+I,J)
00558             END DO
00559          END DO
00560        END IF
00561 *
00562 *      Next Block
00563 *
00564            CUT=CUT+NNB
00565        END DO
00566 *
00567 *        Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H
00568 * 
00569             I=N
00570             DO WHILE ( I .GE. 1 )
00571                IF( IPIV(I) .GT. 0 ) THEN
00572                   IP=IPIV(I)
00573                  IF (I .LT. IP) CALL CHESWAPR( UPLO, N, A, LDA, I ,IP  )
00574                  IF (I .GT. IP) CALL CHESWAPR( UPLO, N, A, LDA, IP ,I )
00575                ELSE
00576                  IP=-IPIV(I)
00577                  IF ( I .LT. IP) CALL CHESWAPR( UPLO, N, A, LDA, I ,IP )
00578                  IF ( I .GT. IP) CALL CHESWAPR( UPLO, N, A, LDA, IP ,I )
00579                  I=I-1
00580                ENDIF
00581                I=I-1
00582             END DO
00583       END IF
00584 *
00585       RETURN
00586 *
00587 *     End of CHETRI2X
00588 *
00589       END
00590 
 All Files Functions