SRC/dsytri.f

   1 *> \brief \b DSYTRI
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSYTRI + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, LDA, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IPIV( * )
  29 *       DOUBLE PRECISION   A( LDA, * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DSYTRI computes the inverse of a real symmetric indefinite matrix
  39 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
  40 *> DSYTRF.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] UPLO
  47 *> \verbatim
  48 *>          UPLO is CHARACTER*1
  49 *>          Specifies whether the details of the factorization are stored
  50 *>          as an upper or lower triangular matrix.
  51 *>          = 'U':  Upper triangular, form is A = U*D*U**T;
  52 *>          = 'L':  Lower triangular, form is A = L*D*L**T.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] N
  56 *> \verbatim
  57 *>          N is INTEGER
  58 *>          The order of the matrix A.  N >= 0.
  59 *> \endverbatim
  60 *>
  61 *> \param[in,out] A
  62 *> \verbatim
  63 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  64 *>          On entry, the block diagonal matrix D and the multipliers
  65 *>          used to obtain the factor U or L as computed by DSYTRF.
  66 *>
  67 *>          On exit, if INFO = 0, the (symmetric) inverse of the original
  68 *>          matrix.  If UPLO = 'U', the upper triangular part of the
  69 *>          inverse is formed and the part of A below the diagonal is not
  70 *>          referenced; if UPLO = 'L' the lower triangular part of the
  71 *>          inverse is formed and the part of A above the diagonal is
  72 *>          not referenced.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] LDA
  76 *> \verbatim
  77 *>          LDA is INTEGER
  78 *>          The leading dimension of the array A.  LDA >= max(1,N).
  79 *> \endverbatim
  80 *>
  81 *> \param[in] IPIV
  82 *> \verbatim
  83 *>          IPIV is INTEGER array, dimension (N)
  84 *>          Details of the interchanges and the block structure of D
  85 *>          as determined by DSYTRF.
  86 *> \endverbatim
  87 *>
  88 *> \param[out] WORK
  89 *> \verbatim
  90 *>          WORK is DOUBLE PRECISION array, dimension (N)
  91 *> \endverbatim
  92 *>
  93 *> \param[out] INFO
  94 *> \verbatim
  95 *>          INFO is INTEGER
  96 *>          = 0: successful exit
  97 *>          < 0: if INFO = -i, the i-th argument had an illegal value
  98 *>          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
  99 *>               inverse could not be computed.
 100 *> \endverbatim
 101 *
 102 *  Authors:
 103 *  ========
 104 *
 105 *> \author Univ. of Tennessee
 106 *> \author Univ. of California Berkeley
 107 *> \author Univ. of Colorado Denver
 108 *> \author NAG Ltd.
 109 *
 110 *> \date November 2011
 111 *
 112 *> \ingroup doubleSYcomputational
 113 *
 114 *  =====================================================================
 115       SUBROUTINE DSYTRI( UPLO, N, A, LDA, IPIV, WORK, INFO )
 116 *
 117 *  -- LAPACK computational routine (version 3.4.0) --
 118 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 119 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 120 *     November 2011
 121 *
 122 *     .. Scalar Arguments ..
 123       CHARACTER          UPLO
 124       INTEGER            INFO, LDA, N
 125 *     ..
 126 *     .. Array Arguments ..
 127       INTEGER            IPIV( * )
 128       DOUBLE PRECISION   A( LDA, * ), WORK( * )
 129 *     ..
 130 *
 131 *  =====================================================================
 132 *
 133 *     .. Parameters ..
 134       DOUBLE PRECISION   ONE, ZERO
 135       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 136 *     ..
 137 *     .. Local Scalars ..
 138       LOGICAL            UPPER
 139       INTEGER            K, KP, KSTEP
 140       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T, TEMP
 141 *     ..
 142 *     .. External Functions ..
 143       LOGICAL            LSAME
 144       DOUBLE PRECISION   DDOT
 145       EXTERNAL           LSAME, DDOT
 146 *     ..
 147 *     .. External Subroutines ..
 148       EXTERNAL           DCOPY, DSWAP, DSYMV, XERBLA
 149 *     ..
 150 *     .. Intrinsic Functions ..
 151       INTRINSIC          ABS, MAX
 152 *     ..
 153 *     .. Executable Statements ..
 154 *
 155 *     Test the input parameters.
 156 *
 157       INFO = 0
 158       UPPER = LSAME( UPLO, 'U' )
 159       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 160          INFO = -1
 161       ELSE IF( N.LT.0 ) THEN
 162          INFO = -2
 163       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 164          INFO = -4
 165       END IF
 166       IF( INFO.NE.0 ) THEN
 167          CALL XERBLA( 'DSYTRI', -INFO )
 168          RETURN
 169       END IF
 170 *
 171 *     Quick return if possible
 172 *
 173       IF( N.EQ.0 )
 174      $   RETURN
 175 *
 176 *     Check that the diagonal matrix D is nonsingular.
 177 *
 178       IF( UPPER ) THEN
 179 *
 180 *        Upper triangular storage: examine D from bottom to top
 181 *
 182          DO 10 INFO = N, 1, -1
 183             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 184      $         RETURN
 185    10    CONTINUE
 186       ELSE
 187 *
 188 *        Lower triangular storage: examine D from top to bottom.
 189 *
 190          DO 20 INFO = 1, N
 191             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 192      $         RETURN
 193    20    CONTINUE
 194       END IF
 195       INFO = 0
 196 *
 197       IF( UPPER ) THEN
 198 *
 199 *        Compute inv(A) from the factorization A = U*D*U**T.
 200 *
 201 *        K is the main loop index, increasing from 1 to N in steps of
 202 *        1 or 2, depending on the size of the diagonal blocks.
 203 *
 204          K = 1
 205    30    CONTINUE
 206 *
 207 *        If K > N, exit from loop.
 208 *
 209          IF( K.GT.N )
 210      $      GO TO 40
 211 *
 212          IF( IPIV( K ).GT.0 ) THEN
 213 *
 214 *           1 x 1 diagonal block
 215 *
 216 *           Invert the diagonal block.
 217 *
 218             A( K, K ) = ONE / A( K, K )
 219 *
 220 *           Compute column K of the inverse.
 221 *
 222             IF( K.GT.1 ) THEN
 223                CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
 224                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
 225      $                     A( 1, K ), 1 )
 226                A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
 227      $                     1 )
 228             END IF
 229             KSTEP = 1
 230          ELSE
 231 *
 232 *           2 x 2 diagonal block
 233 *
 234 *           Invert the diagonal block.
 235 *
 236             T = ABS( A( K, K+1 ) )
 237             AK = A( K, K ) / T
 238             AKP1 = A( K+1, K+1 ) / T
 239             AKKP1 = A( K, K+1 ) / T
 240             D = T*( AK*AKP1-ONE )
 241             A( K, K ) = AKP1 / D
 242             A( K+1, K+1 ) = AK / D
 243             A( K, K+1 ) = -AKKP1 / D
 244 *
 245 *           Compute columns K and K+1 of the inverse.
 246 *
 247             IF( K.GT.1 ) THEN
 248                CALL DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
 249                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
 250      $                     A( 1, K ), 1 )
 251                A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
 252      $                     1 )
 253                A( K, K+1 ) = A( K, K+1 ) -
 254      $                       DDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
 255                CALL DCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
 256                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
 257      $                     A( 1, K+1 ), 1 )
 258                A( K+1, K+1 ) = A( K+1, K+1 ) -
 259      $                         DDOT( K-1, WORK, 1, A( 1, K+1 ), 1 )
 260             END IF
 261             KSTEP = 2
 262          END IF
 263 *
 264          KP = ABS( IPIV( K ) )
 265          IF( KP.NE.K ) THEN
 266 *
 267 *           Interchange rows and columns K and KP in the leading
 268 *           submatrix A(1:k+1,1:k+1)
 269 *
 270             CALL DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
 271             CALL DSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
 272             TEMP = A( K, K )
 273             A( K, K ) = A( KP, KP )
 274             A( KP, KP ) = TEMP
 275             IF( KSTEP.EQ.2 ) THEN
 276                TEMP = A( K, K+1 )
 277                A( K, K+1 ) = A( KP, K+1 )
 278                A( KP, K+1 ) = TEMP
 279             END IF
 280          END IF
 281 *
 282          K = K + KSTEP
 283          GO TO 30
 284    40    CONTINUE
 285 *
 286       ELSE
 287 *
 288 *        Compute inv(A) from the factorization A = L*D*L**T.
 289 *
 290 *        K is the main loop index, increasing from 1 to N in steps of
 291 *        1 or 2, depending on the size of the diagonal blocks.
 292 *
 293          K = N
 294    50    CONTINUE
 295 *
 296 *        If K < 1, exit from loop.
 297 *
 298          IF( K.LT.1 )
 299      $      GO TO 60
 300 *
 301          IF( IPIV( K ).GT.0 ) THEN
 302 *
 303 *           1 x 1 diagonal block
 304 *
 305 *           Invert the diagonal block.
 306 *
 307             A( K, K ) = ONE / A( K, K )
 308 *
 309 *           Compute column K of the inverse.
 310 *
 311             IF( K.LT.N ) THEN
 312                CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
 313                CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
 314      $                     ZERO, A( K+1, K ), 1 )
 315                A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
 316      $                     1 )
 317             END IF
 318             KSTEP = 1
 319          ELSE
 320 *
 321 *           2 x 2 diagonal block
 322 *
 323 *           Invert the diagonal block.
 324 *
 325             T = ABS( A( K, K-1 ) )
 326             AK = A( K-1, K-1 ) / T
 327             AKP1 = A( K, K ) / T
 328             AKKP1 = A( K, K-1 ) / T
 329             D = T*( AK*AKP1-ONE )
 330             A( K-1, K-1 ) = AKP1 / D
 331             A( K, K ) = AK / D
 332             A( K, K-1 ) = -AKKP1 / D
 333 *
 334 *           Compute columns K-1 and K of the inverse.
 335 *
 336             IF( K.LT.N ) THEN
 337                CALL DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
 338                CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
 339      $                     ZERO, A( K+1, K ), 1 )
 340                A( K, K ) = A( K, K ) - DDOT( N-K, WORK, 1, A( K+1, K ),
 341      $                     1 )
 342                A( K, K-1 ) = A( K, K-1 ) -
 343      $                       DDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
 344      $                       1 )
 345                CALL DCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
 346                CALL DSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
 347      $                     ZERO, A( K+1, K-1 ), 1 )
 348                A( K-1, K-1 ) = A( K-1, K-1 ) -
 349      $                         DDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 )
 350             END IF
 351             KSTEP = 2
 352          END IF
 353 *
 354          KP = ABS( IPIV( K ) )
 355          IF( KP.NE.K ) THEN
 356 *
 357 *           Interchange rows and columns K and KP in the trailing
 358 *           submatrix A(k-1:n,k-1:n)
 359 *
 360             IF( KP.LT.N )
 361      $         CALL DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
 362             CALL DSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
 363             TEMP = A( K, K )
 364             A( K, K ) = A( KP, KP )
 365             A( KP, KP ) = TEMP
 366             IF( KSTEP.EQ.2 ) THEN
 367                TEMP = A( K, K-1 )
 368                A( K, K-1 ) = A( KP, K-1 )
 369                A( KP, K-1 ) = TEMP
 370             END IF
 371          END IF
 372 *
 373          K = K - KSTEP
 374          GO TO 50
 375    60    CONTINUE
 376       END IF
 377 *
 378       RETURN
 379 *
 380 *     End of DSYTRI
 381 *
 382       END