LAPACK: SRC/dsptri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSPTRI( UPLO, N, AP, 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, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       DOUBLE PRECISION   AP( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DSPTRI computes the inverse of a real symmetric indefinite matrix
00021 *  A in packed storage using the factorization A = U*D*U**T or
00022 *  A = L*D*L**T computed by DSPTRF.
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 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00037 *          On entry, the block diagonal matrix D and the multipliers
00038 *          used to obtain the factor U or L as computed by DSPTRF,
00039 *          stored as a packed triangular matrix.
00040 *
00041 *          On exit, if INFO = 0, the (symmetric) inverse of the original
00042 *          matrix, stored as a packed triangular matrix. The j-th column
00043 *          of inv(A) is stored in the array AP as follows:
00044 *          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
00045 *          if UPLO = 'L',
00046 *             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
00047 *
00048 *  IPIV    (input) INTEGER array, dimension (N)
00049 *          Details of the interchanges and the block structure of D
00050 *          as determined by DSPTRF.
00051 *
00052 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00053 *
00054 *  INFO    (output) INTEGER
00055 *          = 0: successful exit
00056 *          < 0: if INFO = -i, the i-th argument had an illegal value
00057 *          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
00058 *               inverse could not be computed.
00059 *
00060 *  =====================================================================
00061 *
00062 *     .. Parameters ..
00063       DOUBLE PRECISION   ONE, ZERO
00064       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00065 *     ..
00066 *     .. Local Scalars ..
00067       LOGICAL            UPPER
00068       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
00069       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T, TEMP
00070 *     ..
00071 *     .. External Functions ..
00072       LOGICAL            LSAME
00073       DOUBLE PRECISION   DDOT
00074       EXTERNAL           LSAME, DDOT
00075 *     ..
00076 *     .. External Subroutines ..
00077       EXTERNAL           DCOPY, DSPMV, DSWAP, XERBLA
00078 *     ..
00079 *     .. Intrinsic Functions ..
00080       INTRINSIC          ABS
00081 *     ..
00082 *     .. Executable Statements ..
00083 *
00084 *     Test the input parameters.
00085 *
00086       INFO = 0
00087       UPPER = LSAME( UPLO, 'U' )
00088       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00089          INFO = -1
00090       ELSE IF( N.LT.0 ) THEN
00091          INFO = -2
00092       END IF
00093       IF( INFO.NE.0 ) THEN
00094          CALL XERBLA( 'DSPTRI', -INFO )
00095          RETURN
00096       END IF
00097 *
00098 *     Quick return if possible
00099 *
00100       IF( N.EQ.0 )
00101      $   RETURN
00102 *
00103 *     Check that the diagonal matrix D is nonsingular.
00104 *
00105       IF( UPPER ) THEN
00106 *
00107 *        Upper triangular storage: examine D from bottom to top
00108 *
00109          KP = N*( N+1 ) / 2
00110          DO 10 INFO = N, 1, -1
00111             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
00112      $         RETURN
00113             KP = KP - INFO
00114    10    CONTINUE
00115       ELSE
00116 *
00117 *        Lower triangular storage: examine D from top to bottom.
00118 *
00119          KP = 1
00120          DO 20 INFO = 1, N
00121             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
00122      $         RETURN
00123             KP = KP + N - INFO + 1
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          KC = 1
00137    30    CONTINUE
00138 *
00139 *        If K > N, exit from loop.
00140 *
00141          IF( K.GT.N )
00142      $      GO TO 50
00143 *
00144          KCNEXT = KC + K
00145          IF( IPIV( K ).GT.0 ) THEN
00146 *
00147 *           1 x 1 diagonal block
00148 *
00149 *           Invert the diagonal block.
00150 *
00151             AP( KC+K-1 ) = ONE / AP( KC+K-1 )
00152 *
00153 *           Compute column K of the inverse.
00154 *
00155             IF( K.GT.1 ) THEN
00156                CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
00157                CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
00158      $                     1 )
00159                AP( KC+K-1 ) = AP( KC+K-1 ) -
00160      $                        DDOT( K-1, WORK, 1, AP( KC ), 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( AP( KCNEXT+K-1 ) )
00170             AK = AP( KC+K-1 ) / T
00171             AKP1 = AP( KCNEXT+K ) / T
00172             AKKP1 = AP( KCNEXT+K-1 ) / T
00173             D = T*( AK*AKP1-ONE )
00174             AP( KC+K-1 ) = AKP1 / D
00175             AP( KCNEXT+K ) = AK / D
00176             AP( KCNEXT+K-1 ) = -AKKP1 / D
00177 *
00178 *           Compute columns K and K+1 of the inverse.
00179 *
00180             IF( K.GT.1 ) THEN
00181                CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
00182                CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
00183      $                     1 )
00184                AP( KC+K-1 ) = AP( KC+K-1 ) -
00185      $                        DDOT( K-1, WORK, 1, AP( KC ), 1 )
00186                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
00187      $                            DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
00188      $                            1 )
00189                CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
00190                CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
00191      $                     AP( KCNEXT ), 1 )
00192                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
00193      $                          DDOT( K-1, WORK, 1, AP( KCNEXT ), 1 )
00194             END IF
00195             KSTEP = 2
00196             KCNEXT = KCNEXT + K + 1
00197          END IF
00198 *
00199          KP = ABS( IPIV( K ) )
00200          IF( KP.NE.K ) THEN
00201 *
00202 *           Interchange rows and columns K and KP in the leading
00203 *           submatrix A(1:k+1,1:k+1)
00204 *
00205             KPC = ( KP-1 )*KP / 2 + 1
00206             CALL DSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
00207             KX = KPC + KP - 1
00208             DO 40 J = KP + 1, K - 1
00209                KX = KX + J - 1
00210                TEMP = AP( KC+J-1 )
00211                AP( KC+J-1 ) = AP( KX )
00212                AP( KX ) = TEMP
00213    40       CONTINUE
00214             TEMP = AP( KC+K-1 )
00215             AP( KC+K-1 ) = AP( KPC+KP-1 )
00216             AP( KPC+KP-1 ) = TEMP
00217             IF( KSTEP.EQ.2 ) THEN
00218                TEMP = AP( KC+K+K-1 )
00219                AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
00220                AP( KC+K+KP-1 ) = TEMP
00221             END IF
00222          END IF
00223 *
00224          K = K + KSTEP
00225          KC = KCNEXT
00226          GO TO 30
00227    50    CONTINUE
00228 *
00229       ELSE
00230 *
00231 *        Compute inv(A) from the factorization A = L*D*L'.
00232 *
00233 *        K is the main loop index, increasing from 1 to N in steps of
00234 *        1 or 2, depending on the size of the diagonal blocks.
00235 *
00236          NPP = N*( N+1 ) / 2
00237          K = N
00238          KC = NPP
00239    60    CONTINUE
00240 *
00241 *        If K < 1, exit from loop.
00242 *
00243          IF( K.LT.1 )
00244      $      GO TO 80
00245 *
00246          KCNEXT = KC - ( N-K+2 )
00247          IF( IPIV( K ).GT.0 ) THEN
00248 *
00249 *           1 x 1 diagonal block
00250 *
00251 *           Invert the diagonal block.
00252 *
00253             AP( KC ) = ONE / AP( KC )
00254 *
00255 *           Compute column K of the inverse.
00256 *
00257             IF( K.LT.N ) THEN
00258                CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00259                CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
00260      $                     ZERO, AP( KC+1 ), 1 )
00261                AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
00262             END IF
00263             KSTEP = 1
00264          ELSE
00265 *
00266 *           2 x 2 diagonal block
00267 *
00268 *           Invert the diagonal block.
00269 *
00270             T = ABS( AP( KCNEXT+1 ) )
00271             AK = AP( KCNEXT ) / T
00272             AKP1 = AP( KC ) / T
00273             AKKP1 = AP( KCNEXT+1 ) / T
00274             D = T*( AK*AKP1-ONE )
00275             AP( KCNEXT ) = AKP1 / D
00276             AP( KC ) = AK / D
00277             AP( KCNEXT+1 ) = -AKKP1 / D
00278 *
00279 *           Compute columns K-1 and K of the inverse.
00280 *
00281             IF( K.LT.N ) THEN
00282                CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00283                CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
00284      $                     ZERO, AP( KC+1 ), 1 )
00285                AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
00286                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
00287      $                          DDOT( N-K, AP( KC+1 ), 1,
00288      $                          AP( KCNEXT+2 ), 1 )
00289                CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
00290                CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
00291      $                     ZERO, AP( KCNEXT+2 ), 1 )
00292                AP( KCNEXT ) = AP( KCNEXT ) -
00293      $                        DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
00294             END IF
00295             KSTEP = 2
00296             KCNEXT = KCNEXT - ( N-K+3 )
00297          END IF
00298 *
00299          KP = ABS( IPIV( K ) )
00300          IF( KP.NE.K ) THEN
00301 *
00302 *           Interchange rows and columns K and KP in the trailing
00303 *           submatrix A(k-1:n,k-1:n)
00304 *
00305             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
00306             IF( KP.LT.N )
00307      $         CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
00308             KX = KC + KP - K
00309             DO 70 J = K + 1, KP - 1
00310                KX = KX + N - J + 1
00311                TEMP = AP( KC+J-K )
00312                AP( KC+J-K ) = AP( KX )
00313                AP( KX ) = TEMP
00314    70       CONTINUE
00315             TEMP = AP( KC )
00316             AP( KC ) = AP( KPC )
00317             AP( KPC ) = TEMP
00318             IF( KSTEP.EQ.2 ) THEN
00319                TEMP = AP( KC-N+K-1 )
00320                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
00321                AP( KC-N+KP-1 ) = TEMP
00322             END IF
00323          END IF
00324 *
00325          K = K - KSTEP
00326          KC = KCNEXT
00327          GO TO 60
00328    80    CONTINUE
00329       END IF
00330 *
00331       RETURN
00332 *
00333 *     End of DSPTRI
00334 *
00335       END