LAPACK 3.3.0

chetri.f

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 ),
     $                         1 ) )
00192             END IF
00193             KSTEP = 2
00194          END IF
00195 *
00196          KP = ABS( IPIV( K ) )
00197          IF( KP.NE.K ) THEN
00198 *
00199 *           Interchange rows and columns K and KP in the leading
00200 *           submatrix A(1:k+1,1:k+1)
00201 *
00202             CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
00203             DO 40 J = KP + 1, K - 1
00204                TEMP = CONJG( A( J, K ) )
00205                A( J, K ) = CONJG( A( KP, J ) )
00206                A( KP, J ) = TEMP
00207    40       CONTINUE
00208             A( KP, K ) = CONJG( A( KP, K ) )
00209             TEMP = A( K, K )
00210             A( K, K ) = A( KP, KP )
00211             A( KP, KP ) = TEMP
00212             IF( KSTEP.EQ.2 ) THEN
00213                TEMP = A( K, K+1 )
00214                A( K, K+1 ) = A( KP, K+1 )
00215                A( KP, K+1 ) = TEMP
00216             END IF
00217          END IF
00218 *
00219          K = K + KSTEP
00220          GO TO 30
00221    50    CONTINUE
00222 *
00223       ELSE
00224 *
00225 *        Compute inv(A) from the factorization A = L*D*L'.
00226 *
00227 *        K is the main loop index, increasing from 1 to N in steps of
00228 *        1 or 2, depending on the size of the diagonal blocks.
00229 *
00230          K = N
00231    60    CONTINUE
00232 *
00233 *        If K < 1, exit from loop.
00234 *
00235          IF( K.LT.1 )
00236      $      GO TO 80
00237 *
00238          IF( IPIV( K ).GT.0 ) THEN
00239 *
00240 *           1 x 1 diagonal block
00241 *
00242 *           Invert the diagonal block.
00243 *
00244             A( K, K ) = ONE / REAL( A( K, K ) )
00245 *
00246 *           Compute column K of the inverse.
00247 *
00248             IF( K.LT.N ) THEN
00249                CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00250                CALL CHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
00251      $                     1, ZERO, A( K+1, K ), 1 )
00252                A( K, K ) = A( K, K ) - REAL( CDOTC( N-K, WORK, 1,
     $                     A( K+1, K ), 1 ) )
00253             END IF
00254             KSTEP = 1
00255          ELSE
00256 *
00257 *           2 x 2 diagonal block
00258 *
00259 *           Invert the diagonal block.
00260 *
00261             T = ABS( A( K, K-1 ) )
00262             AK = REAL( A( K-1, K-1 ) ) / T
00263             AKP1 = REAL( A( K, K ) ) / T
00264             AKKP1 = A( K, K-1 ) / T
00265             D = T*( AK*AKP1-ONE )
00266             A( K-1, K-1 ) = AKP1 / D
00267             A( K, K ) = AK / D
00268             A( K, K-1 ) = -AKKP1 / D
00269 *
00270 *           Compute columns K-1 and K of the inverse.
00271 *
00272             IF( K.LT.N ) THEN
00273                CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00274                CALL CHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
00275      $                     1, ZERO, A( K+1, K ), 1 )
00276                A( K, K ) = A( K, K ) - REAL( CDOTC( N-K, WORK, 1,
     $                     A( K+1, K ), 1 ) )
00277                A( K, K-1 ) = A( K, K-1 ) -
00278      $                       CDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
00279      $                       1 )
00280                CALL CCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
00281                CALL CHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
00282      $                     1, ZERO, A( K+1, K-1 ), 1 )
00283                A( K-1, K-1 ) = A( K-1, K-1 ) -
00284      $                         REAL( CDOTC( N-K, WORK, 1, A( K+1, K-1 ),
     $                         1 ) )
00285             END IF
00286             KSTEP = 2
00287          END IF
00288 *
00289          KP = ABS( IPIV( K ) )
00290          IF( KP.NE.K ) THEN
00291 *
00292 *           Interchange rows and columns K and KP in the trailing
00293 *           submatrix A(k-1:n,k-1:n)
00294 *
00295             IF( KP.LT.N )
00296      $         CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
00297             DO 70 J = K + 1, KP - 1
00298                TEMP = CONJG( A( J, K ) )
00299                A( J, K ) = CONJG( A( KP, J ) )
00300                A( KP, J ) = TEMP
00301    70       CONTINUE
00302             A( KP, K ) = CONJG( A( KP, K ) )
00303             TEMP = A( K, K )
00304             A( K, K ) = A( KP, KP )
00305             A( KP, KP ) = TEMP
00306             IF( KSTEP.EQ.2 ) THEN
00307                TEMP = A( K, K-1 )
00308                A( K, K-1 ) = A( KP, K-1 )
00309                A( KP, K-1 ) = TEMP
00310             END IF
00311          END IF
00312 *
00313          K = K - KSTEP
00314          GO TO 60
00315    80    CONTINUE
00316       END IF
00317 *
00318       RETURN
00319 *
00320 *     End of CHETRI
00321 *
00322       END
00323 
 All Files Functions