SRC/csptri.f

   1 *> \brief \b CSPTRI
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CSPTRI + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csptri.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csptri.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csptri.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CSPTRI( 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 *       COMPLEX            AP( * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CSPTRI computes the inverse of a complex 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 CSPTRF.
  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 COMPLEX 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 CSPTRF,
  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 CSPTRF.
  81 *> \endverbatim
  82 *>
  83 *> \param[out] WORK
  84 *> \verbatim
  85 *>          WORK is COMPLEX 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 complexOTHERcomputational
 108 *
 109 *  =====================================================================
 110       SUBROUTINE CSPTRI( 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       COMPLEX            AP( * ), WORK( * )
 124 *     ..
 125 *
 126 *  =====================================================================
 127 *
 128 *     .. Parameters ..
 129       COMPLEX            ONE, ZERO
 130       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
 131      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
 132 *     ..
 133 *     .. Local Scalars ..
 134       LOGICAL            UPPER
 135       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
 136       COMPLEX            AK, AKKP1, AKP1, D, T, TEMP
 137 *     ..
 138 *     .. External Functions ..
 139       LOGICAL            LSAME
 140       COMPLEX            CDOTU
 141       EXTERNAL           LSAME, CDOTU
 142 *     ..
 143 *     .. External Subroutines ..
 144       EXTERNAL           CCOPY, CSPMV, CSWAP, XERBLA
 145 *     ..
 146 *     .. Intrinsic Functions ..
 147       INTRINSIC          ABS
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Test the input parameters.
 152 *
 153       INFO = 0
 154       UPPER = LSAME( UPLO, 'U' )
 155       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 156          INFO = -1
 157       ELSE IF( N.LT.0 ) THEN
 158          INFO = -2
 159       END IF
 160       IF( INFO.NE.0 ) THEN
 161          CALL XERBLA( 'CSPTRI', -INFO )
 162          RETURN
 163       END IF
 164 *
 165 *     Quick return if possible
 166 *
 167       IF( N.EQ.0 )
 168      $   RETURN
 169 *
 170 *     Check that the diagonal matrix D is nonsingular.
 171 *
 172       IF( UPPER ) THEN
 173 *
 174 *        Upper triangular storage: examine D from bottom to top
 175 *
 176          KP = N*( N+1 ) / 2
 177          DO 10 INFO = N, 1, -1
 178             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
 179      $         RETURN
 180             KP = KP - INFO
 181    10    CONTINUE
 182       ELSE
 183 *
 184 *        Lower triangular storage: examine D from top to bottom.
 185 *
 186          KP = 1
 187          DO 20 INFO = 1, N
 188             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
 189      $         RETURN
 190             KP = KP + N - INFO + 1
 191    20    CONTINUE
 192       END IF
 193       INFO = 0
 194 *
 195       IF( UPPER ) THEN
 196 *
 197 *        Compute inv(A) from the factorization A = U*D*U**T.
 198 *
 199 *        K is the main loop index, increasing from 1 to N in steps of
 200 *        1 or 2, depending on the size of the diagonal blocks.
 201 *
 202          K = 1
 203          KC = 1
 204    30    CONTINUE
 205 *
 206 *        If K > N, exit from loop.
 207 *
 208          IF( K.GT.N )
 209      $      GO TO 50
 210 *
 211          KCNEXT = KC + K
 212          IF( IPIV( K ).GT.0 ) THEN
 213 *
 214 *           1 x 1 diagonal block
 215 *
 216 *           Invert the diagonal block.
 217 *
 218             AP( KC+K-1 ) = ONE / AP( KC+K-1 )
 219 *
 220 *           Compute column K of the inverse.
 221 *
 222             IF( K.GT.1 ) THEN
 223                CALL CCOPY( K-1, AP( KC ), 1, WORK, 1 )
 224                CALL CSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
 225      $                     1 )
 226                AP( KC+K-1 ) = AP( KC+K-1 ) -
 227      $                        CDOTU( K-1, WORK, 1, AP( KC ), 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 = AP( KCNEXT+K-1 )
 237             AK = AP( KC+K-1 ) / T
 238             AKP1 = AP( KCNEXT+K ) / T
 239             AKKP1 = AP( KCNEXT+K-1 ) / T
 240             D = T*( AK*AKP1-ONE )
 241             AP( KC+K-1 ) = AKP1 / D
 242             AP( KCNEXT+K ) = AK / D
 243             AP( KCNEXT+K-1 ) = -AKKP1 / D
 244 *
 245 *           Compute columns K and K+1 of the inverse.
 246 *
 247             IF( K.GT.1 ) THEN
 248                CALL CCOPY( K-1, AP( KC ), 1, WORK, 1 )
 249                CALL CSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
 250      $                     1 )
 251                AP( KC+K-1 ) = AP( KC+K-1 ) -
 252      $                        CDOTU( K-1, WORK, 1, AP( KC ), 1 )
 253                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
 254      $                            CDOTU( K-1, AP( KC ), 1, AP( KCNEXT ),
 255      $                            1 )
 256                CALL CCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
 257                CALL CSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
 258      $                     AP( KCNEXT ), 1 )
 259                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
 260      $                          CDOTU( K-1, WORK, 1, AP( KCNEXT ), 1 )
 261             END IF
 262             KSTEP = 2
 263             KCNEXT = KCNEXT + K + 1
 264          END IF
 265 *
 266          KP = ABS( IPIV( K ) )
 267          IF( KP.NE.K ) THEN
 268 *
 269 *           Interchange rows and columns K and KP in the leading
 270 *           submatrix A(1:k+1,1:k+1)
 271 *
 272             KPC = ( KP-1 )*KP / 2 + 1
 273             CALL CSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
 274             KX = KPC + KP - 1
 275             DO 40 J = KP + 1, K - 1
 276                KX = KX + J - 1
 277                TEMP = AP( KC+J-1 )
 278                AP( KC+J-1 ) = AP( KX )
 279                AP( KX ) = TEMP
 280    40       CONTINUE
 281             TEMP = AP( KC+K-1 )
 282             AP( KC+K-1 ) = AP( KPC+KP-1 )
 283             AP( KPC+KP-1 ) = TEMP
 284             IF( KSTEP.EQ.2 ) THEN
 285                TEMP = AP( KC+K+K-1 )
 286                AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
 287                AP( KC+K+KP-1 ) = TEMP
 288             END IF
 289          END IF
 290 *
 291          K = K + KSTEP
 292          KC = KCNEXT
 293          GO TO 30
 294    50    CONTINUE
 295 *
 296       ELSE
 297 *
 298 *        Compute inv(A) from the factorization A = L*D*L**T.
 299 *
 300 *        K is the main loop index, increasing from 1 to N in steps of
 301 *        1 or 2, depending on the size of the diagonal blocks.
 302 *
 303          NPP = N*( N+1 ) / 2
 304          K = N
 305          KC = NPP
 306    60    CONTINUE
 307 *
 308 *        If K < 1, exit from loop.
 309 *
 310          IF( K.LT.1 )
 311      $      GO TO 80
 312 *
 313          KCNEXT = KC - ( N-K+2 )
 314          IF( IPIV( K ).GT.0 ) THEN
 315 *
 316 *           1 x 1 diagonal block
 317 *
 318 *           Invert the diagonal block.
 319 *
 320             AP( KC ) = ONE / AP( KC )
 321 *
 322 *           Compute column K of the inverse.
 323 *
 324             IF( K.LT.N ) THEN
 325                CALL CCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
 326                CALL CSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
 327      $                     ZERO, AP( KC+1 ), 1 )
 328                AP( KC ) = AP( KC ) - CDOTU( N-K, WORK, 1, AP( KC+1 ),
 329      $                    1 )
 330             END IF
 331             KSTEP = 1
 332          ELSE
 333 *
 334 *           2 x 2 diagonal block
 335 *
 336 *           Invert the diagonal block.
 337 *
 338             T = AP( KCNEXT+1 )
 339             AK = AP( KCNEXT ) / T
 340             AKP1 = AP( KC ) / T
 341             AKKP1 = AP( KCNEXT+1 ) / T
 342             D = T*( AK*AKP1-ONE )
 343             AP( KCNEXT ) = AKP1 / D
 344             AP( KC ) = AK / D
 345             AP( KCNEXT+1 ) = -AKKP1 / D
 346 *
 347 *           Compute columns K-1 and K of the inverse.
 348 *
 349             IF( K.LT.N ) THEN
 350                CALL CCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
 351                CALL CSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
 352      $                     ZERO, AP( KC+1 ), 1 )
 353                AP( KC ) = AP( KC ) - CDOTU( N-K, WORK, 1, AP( KC+1 ),
 354      $                    1 )
 355                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
 356      $                          CDOTU( N-K, AP( KC+1 ), 1,
 357      $                          AP( KCNEXT+2 ), 1 )
 358                CALL CCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
 359                CALL CSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
 360      $                     ZERO, AP( KCNEXT+2 ), 1 )
 361                AP( KCNEXT ) = AP( KCNEXT ) -
 362      $                        CDOTU( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
 363             END IF
 364             KSTEP = 2
 365             KCNEXT = KCNEXT - ( N-K+3 )
 366          END IF
 367 *
 368          KP = ABS( IPIV( K ) )
 369          IF( KP.NE.K ) THEN
 370 *
 371 *           Interchange rows and columns K and KP in the trailing
 372 *           submatrix A(k-1:n,k-1:n)
 373 *
 374             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
 375             IF( KP.LT.N )
 376      $         CALL CSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
 377             KX = KC + KP - K
 378             DO 70 J = K + 1, KP - 1
 379                KX = KX + N - J + 1
 380                TEMP = AP( KC+J-K )
 381                AP( KC+J-K ) = AP( KX )
 382                AP( KX ) = TEMP
 383    70       CONTINUE
 384             TEMP = AP( KC )
 385             AP( KC ) = AP( KPC )
 386             AP( KPC ) = TEMP
 387             IF( KSTEP.EQ.2 ) THEN
 388                TEMP = AP( KC-N+K-1 )
 389                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
 390                AP( KC-N+KP-1 ) = TEMP
 391             END IF
 392          END IF
 393 *
 394          K = K - KSTEP
 395          KC = KCNEXT
 396          GO TO 60
 397    80    CONTINUE
 398       END IF
 399 *
 400       RETURN
 401 *
 402 *     End of CSPTRI
 403 *
 404       END