LAPACK 3.3.1
Linear Algebra PACKage

chetri2x.f

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