LAPACK: SRC/chetri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHETRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, LDA, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       COMPLEX            A( LDA, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CHETRI computes the inverse of a complex Hermitian indefinite matrix
00021 *  A using the factorization A = U*D*U**H or A = L*D*L**H computed by
00022 *  CHETRF.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  UPLO    (input) CHARACTER*1
00028 *          Specifies whether the details of the factorization are stored
00029 *          as an upper or lower triangular matrix.
00030 *          = 'U':  Upper triangular, form is A = U*D*U**H;
00031 *          = 'L':  Lower triangular, form is A = L*D*L**H.
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00037 *          On entry, the block diagonal matrix D and the multipliers
00038 *          used to obtain the factor U or L as computed by CHETRF.
00039 *
00040 *          On exit, if INFO = 0, the (Hermitian) inverse of the original
00041 *          matrix.  If UPLO = 'U', the upper triangular part of the
00042 *          inverse is formed and the part of A below the diagonal is not
00043 *          referenced; if UPLO = 'L' the lower triangular part of the
00044 *          inverse is formed and the part of A above the diagonal is
00045 *          not referenced.
00046 *
00047 *  LDA     (input) INTEGER
00048 *          The leading dimension of the array A.  LDA >= max(1,N).
00049 *
00050 *  IPIV    (input) INTEGER array, dimension (N)
00051 *          Details of the interchanges and the block structure of D
00052 *          as determined by CHETRF.
00053 *
00054 *  WORK    (workspace) COMPLEX array, dimension (N)
00055 *
00056 *  INFO    (output) INTEGER
00057 *          = 0: successful exit
00058 *          < 0: if INFO = -i, the i-th argument had an illegal value
00059 *          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
00060 *               inverse could not be computed.
00061 *
00062 *  =====================================================================
00063 *
00064 *     .. Parameters ..
00065       REAL               ONE
00066       COMPLEX            CONE, ZERO
00067       PARAMETER          ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ),
00068      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
00069 *     ..
00070 *     .. Local Scalars ..
00071       LOGICAL            UPPER
00072       INTEGER            J, K, KP, KSTEP
00073       REAL               AK, AKP1, D, T
00074       COMPLEX            AKKP1, TEMP
00075 *     ..
00076 *     .. External Functions ..
00077       LOGICAL            LSAME
00078       COMPLEX            CDOTC
00079       EXTERNAL           LSAME, CDOTC
00080 *     ..
00081 *     .. External Subroutines ..
00082       EXTERNAL           CCOPY, CHEMV, CSWAP, XERBLA
00083 *     ..
00084 *     .. Intrinsic Functions ..
00085       INTRINSIC          ABS, CONJG, MAX, REAL
00086 *     ..
00087 *     .. Executable Statements ..
00088 *
00089 *     Test the input parameters.
00090 *
00091       INFO = 0
00092       UPPER = LSAME( UPLO, 'U' )
00093       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00094          INFO = -1
00095       ELSE IF( N.LT.0 ) THEN
00096          INFO = -2
00097       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00098          INFO = -4
00099       END IF
00100       IF( INFO.NE.0 ) THEN
00101          CALL XERBLA( 'CHETRI', -INFO )
00102          RETURN
00103       END IF
00104 *
00105 *     Quick return if possible
00106 *
00107       IF( N.EQ.0 )
00108      $   RETURN
00109 *
00110 *     Check that the diagonal matrix D is nonsingular.
00111 *
00112       IF( UPPER ) THEN
00113 *
00114 *        Upper triangular storage: examine D from bottom to top
00115 *
00116          DO 10 INFO = N, 1, -1
00117             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00118      $         RETURN
00119    10    CONTINUE
00120       ELSE
00121 *
00122 *        Lower triangular storage: examine D from top to bottom.
00123 *
00124          DO 20 INFO = 1, N
00125             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00126      $         RETURN
00127    20    CONTINUE
00128       END IF
00129       INFO = 0
00130 *
00131       IF( UPPER ) THEN
00132 *
00133 *        Compute inv(A) from the factorization A = U*D*U'.
00134 *
00135 *        K is the main loop index, increasing from 1 to N in steps of
00136 *        1 or 2, depending on the size of the diagonal blocks.
00137 *
00138          K = 1
00139    30    CONTINUE
00140 *
00141 *        If K > N, exit from loop.
00142 *
00143          IF( K.GT.N )
00144      $      GO TO 50
00145 *
00146          IF( IPIV( K ).GT.0 ) THEN
00147 *
00148 *           1 x 1 diagonal block
00149 *
00150 *           Invert the diagonal block.
00151 *
00152             A( K, K ) = ONE / REAL( A( K, K ) )
00153 *
00154 *           Compute column K of the inverse.
00155 *
00156             IF( K.GT.1 ) THEN
00157                CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 )
00158                CALL CHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
00159      $                     A( 1, K ), 1 )
00160                A( K, K ) = A( K, K ) - REAL( CDOTC( K-1, WORK, 1, A( 1,
     $                     K ), 1 ) )
00161             END IF
00162             KSTEP = 1
00163          ELSE
00164 *
00165 *           2 x 2 diagonal block
00166 *
00167 *           Invert the diagonal block.
00168 *
00169             T = ABS( A( K, K+1 ) )
00170             AK = REAL( A( K, K ) ) / T
00171             AKP1 = REAL( A( K+1, K+1 ) ) / T
00172             AKKP1 = A( K, K+1 ) / T
00173             D = T*( AK*AKP1-ONE )
00174             A( K, K ) = AKP1 / D
00175             A( K+1, K+1 ) = AK / D
00176             A( K, K+1 ) = -AKKP1 / D
00177 *
00178 *           Compute columns K and K+1 of the inverse.
00179 *
00180             IF( K.GT.1 ) THEN
00181                CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 )
00182                CALL CHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
00183      $                     A( 1, K ), 1 )
00184                A( K, K ) = A( K, K ) - REAL( CDOTC( K-1, WORK, 1, A( 1,
     $                     K ), 1 ) )
00185                A( K, K+1 ) = A( K, K+1 ) -
00186      $                       CDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
00187                CALL CCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
00188                CALL CHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, ZERO,
00189      $                     A( 1, K+1 ), 1 )
00190                A( K+1, K+1 ) = A( K+1, K+1 ) -
00191      $                         REAL( CDOTC( K-1, WORK, 1, A( 1, K+1 ),
00192      $                         1 ) )
00193             END IF
00194             KSTEP = 2
00195          END IF
00196 *
00197          KP = ABS( IPIV( K ) )
00198          IF( KP.NE.K ) THEN
00199 *
00200 *           Interchange rows and columns K and KP in the leading
00201 *           submatrix A(1:k+1,1:k+1)
00202 *
00203             CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
00204             DO 40 J = KP + 1, K - 1
00205                TEMP = CONJG( A( J, K ) )
00206                A( J, K ) = CONJG( A( KP, J ) )
00207                A( KP, J ) = TEMP
00208    40       CONTINUE
00209             A( KP, K ) = CONJG( A( KP, K ) )
00210             TEMP = A( K, K )
00211             A( K, K ) = A( KP, KP )
00212             A( KP, KP ) = TEMP
00213             IF( KSTEP.EQ.2 ) THEN
00214                TEMP = A( K, K+1 )
00215                A( K, K+1 ) = A( KP, K+1 )
00216                A( KP, K+1 ) = TEMP
00217             END IF
00218          END IF
00219 *
00220          K = K + KSTEP
00221          GO TO 30
00222    50    CONTINUE
00223 *
00224       ELSE
00225 *
00226 *        Compute inv(A) from the factorization A = L*D*L'.
00227 *
00228 *        K is the main loop index, increasing from 1 to N in steps of
00229 *        1 or 2, depending on the size of the diagonal blocks.
00230 *
00231          K = N
00232    60    CONTINUE
00233 *
00234 *        If K < 1, exit from loop.
00235 *
00236          IF( K.LT.1 )
00237      $      GO TO 80
00238 *
00239          IF( IPIV( K ).GT.0 ) THEN
00240 *
00241 *           1 x 1 diagonal block
00242 *
00243 *           Invert the diagonal block.
00244 *
00245             A( K, K ) = ONE / REAL( A( K, K ) )
00246 *
00247 *           Compute column K of the inverse.
00248 *
00249             IF( K.LT.N ) THEN
00250                CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00251                CALL CHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
00252      $                     1, ZERO, A( K+1, K ), 1 )
00253                A( K, K ) = A( K, K ) - REAL( CDOTC( N-K, WORK, 1,
     $                     A( K+1, K ), 1 ) )
00254             END IF
00255             KSTEP = 1
00256          ELSE
00257 *
00258 *           2 x 2 diagonal block
00259 *
00260 *           Invert the diagonal block.
00261 *
00262             T = ABS( A( K, K-1 ) )
00263             AK = REAL( A( K-1, K-1 ) ) / T
00264             AKP1 = REAL( A( K, K ) ) / T
00265             AKKP1 = A( K, K-1 ) / T
00266             D = T*( AK*AKP1-ONE )
00267             A( K-1, K-1 ) = AKP1 / D
00268             A( K, K ) = AK / D
00269             A( K, K-1 ) = -AKKP1 / D
00270 *
00271 *           Compute columns K-1 and K of the inverse.
00272 *
00273             IF( K.LT.N ) THEN
00274                CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00275                CALL CHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
00276      $                     1, ZERO, A( K+1, K ), 1 )
00277                A( K, K ) = A( K, K ) - REAL( CDOTC( N-K, WORK, 1,
     $                     A( K+1, K ), 1 ) )
00278                A( K, K-1 ) = A( K, K-1 ) -
00279      $                       CDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
00280      $                       1 )
00281                CALL CCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
00282                CALL CHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
00283      $                     1, ZERO, A( K+1, K-1 ), 1 )
00284                A( K-1, K-1 ) = A( K-1, K-1 ) -
00285      $                         REAL( CDOTC( N-K, WORK, 1, A( K+1, K-1 ),
00286      $                         1 ) )
00287             END IF
00288             KSTEP = 2
00289          END IF
00290 *
00291          KP = ABS( IPIV( K ) )
00292          IF( KP.NE.K ) THEN
00293 *
00294 *           Interchange rows and columns K and KP in the trailing
00295 *           submatrix A(k-1:n,k-1:n)
00296 *
00297             IF( KP.LT.N )
00298      $         CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
00299             DO 70 J = K + 1, KP - 1
00300                TEMP = CONJG( A( J, K ) )
00301                A( J, K ) = CONJG( A( KP, J ) )
00302                A( KP, J ) = TEMP
00303    70       CONTINUE
00304             A( KP, K ) = CONJG( A( KP, K ) )
00305             TEMP = A( K, K )
00306             A( K, K ) = A( KP, KP )
00307             A( KP, KP ) = TEMP
00308             IF( KSTEP.EQ.2 ) THEN
00309                TEMP = A( K, K-1 )
00310                A( K, K-1 ) = A( KP, K-1 )
00311                A( KP, K-1 ) = TEMP
00312             END IF
00313          END IF
00314 *
00315          K = K - KSTEP
00316          GO TO 60
00317    80    CONTINUE
00318       END IF
00319 *
00320       RETURN
00321 *
00322 *     End of CHETRI
00323 *
00324       END
00325