LAPACK: SRC/dsytri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSYTRI( 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       DOUBLE PRECISION   A( LDA, * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DSYTRI computes the inverse of a real symmetric indefinite matrix
00021 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
00022 *  DSYTRF.
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**T;
00031 *          = 'L':  Lower triangular, form is A = L*D*L**T.
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  A       (input/output) DOUBLE PRECISION 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 DSYTRF.
00039 *
00040 *          On exit, if INFO = 0, the (symmetric) 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 DSYTRF.
00053 *
00054 *  WORK    (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ONE, ZERO
00066       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00067 *     ..
00068 *     .. Local Scalars ..
00069       LOGICAL            UPPER
00070       INTEGER            K, KP, KSTEP
00071       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T, TEMP
00072 *     ..
00073 *     .. External Functions ..
00074       LOGICAL            LSAME
00075       DOUBLE PRECISION   DDOT
00076       EXTERNAL           LSAME, DDOT
00077 *     ..
00078 *     .. External Subroutines ..
00079       EXTERNAL           DCOPY, DSWAP, DSYMV, XERBLA
00080 *     ..
00081 *     .. Intrinsic Functions ..
00082       INTRINSIC          ABS, MAX
00083 *     ..
00084 *     .. Executable Statements ..
00085 *
00086 *     Test the input parameters.
00087 *
00088       INFO = 0
00089       UPPER = LSAME( UPLO, 'U' )
00090       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00091          INFO = -1
00092       ELSE IF( N.LT.0 ) THEN
00093          INFO = -2
00094       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00095          INFO = -4
00096       END IF
00097       IF( INFO.NE.0 ) THEN
00098          CALL XERBLA( 'DSYTRI', -INFO )
00099          RETURN
00100       END IF
00101 *
00102 *     Quick return if possible
00103 *
00104       IF( N.EQ.0 )
00105      $   RETURN
00106 *
00107 *     Check that the diagonal matrix D is nonsingular.
00108 *
00109       IF( UPPER ) THEN
00110 *
00111 *        Upper triangular storage: examine D from bottom to top
00112 *
00113          DO 10 INFO = N, 1, -1
00114             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00115      $         RETURN
00116    10    CONTINUE
00117       ELSE
00118 *
00119 *        Lower triangular storage: examine D from top to bottom.
00120 *
00121          DO 20 INFO = 1, N
00122             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00123      $         RETURN
00124    20    CONTINUE
00125       END IF
00126       INFO = 0
00127 *
00128       IF( UPPER ) THEN
00129 *
00130 *        Compute inv(A) from the factorization A = U*D*U'.
00131 *
00132 *        K is the main loop index, increasing from 1 to N in steps of
00133 *        1 or 2, depending on the size of the diagonal blocks.
00134 *
00135          K = 1
00136    30    CONTINUE
00137 *
00138 *        If K > N, exit from loop.
00139 *
00140          IF( K.GT.N )
00141      $      GO TO 40
00142 *
00143          IF( IPIV( K ).GT.0 ) THEN
00144 *
00145 *           1 x 1 diagonal block
00146 *
00147 *           Invert the diagonal block.
00148 *
00149             A( K, K ) = ONE / A( K, K )
00150 *
00151 *           Compute column K of the inverse.
00152 *
00153             IF( K.GT.1 ) THEN
00154                CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
00155                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
00156      $                     A( 1, K ), 1 )
00157                A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
00158      $                     1 )
00159             END IF
00160             KSTEP = 1
00161          ELSE
00162 *
00163 *           2 x 2 diagonal block
00164 *
00165 *           Invert the diagonal block.
00166 *
00167             T = ABS( A( K, K+1 ) )
00168             AK = A( K, K ) / T
00169             AKP1 = A( K+1, K+1 ) / T
00170             AKKP1 = A( K, K+1 ) / T
00171             D = T*( AK*AKP1-ONE )
00172             A( K, K ) = AKP1 / D
00173             A( K+1, K+1 ) = AK / D
00174             A( K, K+1 ) = -AKKP1 / D
00175 *
00176 *           Compute columns K and K+1 of the inverse.
00177 *
00178             IF( K.GT.1 ) THEN
00179                CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
00180                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
00181      $                     A( 1, K ), 1 )
00182                A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
00183      $                     1 )
00184                A( K, K+1 ) = A( K, K+1 ) -
00185      $                       DDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
00186                CALL DCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
00187                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
00188      $                     A( 1, K+1 ), 1 )
00189                A( K+1, K+1 ) = A( K+1, K+1 ) -
00190      $                         DDOT( K-1, WORK, 1, A( 1, K+1 ), 1 )
00191             END IF
00192             KSTEP = 2
00193          END IF
00194 *
00195          KP = ABS( IPIV( K ) )
00196          IF( KP.NE.K ) THEN
00197 *
00198 *           Interchange rows and columns K and KP in the leading
00199 *           submatrix A(1:k+1,1:k+1)
00200 *
00201             CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
00202             CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
00203             TEMP = A( K, K )
00204             A( K, K ) = A( KP, KP )
00205             A( KP, KP ) = TEMP
00206             IF( KSTEP.EQ.2 ) THEN
00207                TEMP = A( K, K+1 )
00208                A( K, K+1 ) = A( KP, K+1 )
00209                A( KP, K+1 ) = TEMP
00210             END IF
00211          END IF
00212 *
00213          K = K + KSTEP
00214          GO TO 30
00215    40    CONTINUE
00216 *
00217       ELSE
00218 *
00219 *        Compute inv(A) from the factorization A = L*D*L'.
00220 *
00221 *        K is the main loop index, increasing from 1 to N in steps of
00222 *        1 or 2, depending on the size of the diagonal blocks.
00223 *
00224          K = N
00225    50    CONTINUE
00226 *
00227 *        If K < 1, exit from loop.
00228 *
00229          IF( K.LT.1 )
00230      $      GO TO 60
00231 *
00232          IF( IPIV( K ).GT.0 ) THEN
00233 *
00234 *           1 x 1 diagonal block
00235 *
00236 *           Invert the diagonal block.
00237 *
00238             A( K, K ) = ONE / A( K, K )
00239 *
00240 *           Compute column K of the inverse.
00241 *
00242             IF( K.LT.N ) THEN
00243                CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00244                CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
00245      $                     ZERO, A( K+1, K ), 1 )
00246                A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
00247      $                     1 )
00248             END IF
00249             KSTEP = 1
00250          ELSE
00251 *
00252 *           2 x 2 diagonal block
00253 *
00254 *           Invert the diagonal block.
00255 *
00256             T = ABS( A( K, K-1 ) )
00257             AK = A( K-1, K-1 ) / T
00258             AKP1 = A( K, K ) / T
00259             AKKP1 = A( K, K-1 ) / T
00260             D = T*( AK*AKP1-ONE )
00261             A( K-1, K-1 ) = AKP1 / D
00262             A( K, K ) = AK / D
00263             A( K, K-1 ) = -AKKP1 / D
00264 *
00265 *           Compute columns K-1 and K of the inverse.
00266 *
00267             IF( K.LT.N ) THEN
00268                CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
00269                CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
00270      $                     ZERO, A( K+1, K ), 1 )
00271                A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
00272      $                     1 )
00273                A( K, K-1 ) = A( K, K-1 ) -
00274      $                       DDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
00275      $                       1 )
00276                CALL DCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
00277                CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
00278      $                     ZERO, A( K+1, K-1 ), 1 )
00279                A( K-1, K-1 ) = A( K-1, K-1 ) -
00280      $                         DDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 )
00281             END IF
00282             KSTEP = 2
00283          END IF
00284 *
00285          KP = ABS( IPIV( K ) )
00286          IF( KP.NE.K ) THEN
00287 *
00288 *           Interchange rows and columns K and KP in the trailing
00289 *           submatrix A(k-1:n,k-1:n)
00290 *
00291             IF( KP.LT.N )
00292      $         CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
00293             CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
00294             TEMP = A( K, K )
00295             A( K, K ) = A( KP, KP )
00296             A( KP, KP ) = TEMP
00297             IF( KSTEP.EQ.2 ) THEN
00298                TEMP = A( K, K-1 )
00299                A( K, K-1 ) = A( KP, K-1 )
00300                A( KP, K-1 ) = TEMP
00301             END IF
00302          END IF
00303 *
00304          K = K - KSTEP
00305          GO TO 50
00306    60    CONTINUE
00307       END IF
00308 *
00309       RETURN
00310 *
00311 *     End of DSYTRI
00312 *
00313       END