SRC/dsptri.f

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