SRC/ssytri_rook.f

   1 *> \brief \b SSYTRI_ROOK
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSYTRI_ROOK + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri_rook.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri_rook.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri_rook.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSYTRI_ROOK( 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 *       REAL               A( LDA, * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SSYTRI_ROOK computes the inverse of a real symmetric
  39 *> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
  40 *> computed by SSYTRF_ROOK.
  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 REAL 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 SSYTRF_ROOK.
  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 SSYTRF_ROOK.
  86 *> \endverbatim
  87 *>
  88 *> \param[out] WORK
  89 *> \verbatim
  90 *>          WORK is REAL 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 April 2012
 111 *
 112 *> \ingroup realSYcomputational
 113 *
 114 *> \par Contributors:
 115 *  ==================
 116 *>
 117 *> \verbatim
 118 *>
 119 *>   April 2012, Igor Kozachenko,
 120 *>                  Computer Science Division,
 121 *>                  University of California, Berkeley
 122 *>
 123 *>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
 124 *>                  School of Mathematics,
 125 *>                  University of Manchester
 126 *>
 127 *> \endverbatim
 128 *
 129 *  =====================================================================
 130       SUBROUTINE SSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
 131 *
 132 *  -- LAPACK computational routine (version 3.4.1) --
 133 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 134 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 135 *     April 2012
 136 *
 137 *     .. Scalar Arguments ..
 138       CHARACTER          UPLO
 139       INTEGER            INFO, LDA, N
 140 *     ..
 141 *     .. Array Arguments ..
 142       INTEGER            IPIV( * )
 143       REAL               A( LDA, * ), WORK( * )
 144 *     ..
 145 *
 146 *  =====================================================================
 147 *
 148 *     .. Parameters ..
 149       REAL               ONE, ZERO
 150       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 151 *     ..
 152 *     .. Local Scalars ..
 153       LOGICAL            UPPER
 154       INTEGER            K, KP, KSTEP
 155       REAL               AK, AKKP1, AKP1, D, T, TEMP
 156 *     ..
 157 *     .. External Functions ..
 158       LOGICAL            LSAME
 159       REAL               SDOT
 160       EXTERNAL           LSAME, SDOT
 161 *     ..
 162 *     .. External Subroutines ..
 163       EXTERNAL           SCOPY, SSWAP, SSYMV, XERBLA
 164 *     ..
 165 *     .. Intrinsic Functions ..
 166       INTRINSIC          ABS, MAX
 167 *     ..
 168 *     .. Executable Statements ..
 169 *
 170 *     Test the input parameters.
 171 *
 172       INFO = 0
 173       UPPER = LSAME( UPLO, 'U' )
 174       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 175          INFO = -1
 176       ELSE IF( N.LT.0 ) THEN
 177          INFO = -2
 178       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 179          INFO = -4
 180       END IF
 181       IF( INFO.NE.0 ) THEN
 182          CALL XERBLA( 'SSYTRI_ROOK', -INFO )
 183          RETURN
 184       END IF
 185 *
 186 *     Quick return if possible
 187 *
 188       IF( N.EQ.0 )
 189      $   RETURN
 190 *
 191 *     Check that the diagonal matrix D is nonsingular.
 192 *
 193       IF( UPPER ) THEN
 194 *
 195 *        Upper triangular storage: examine D from bottom to top
 196 *
 197          DO 10 INFO = N, 1, -1
 198             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 199      $         RETURN
 200    10    CONTINUE
 201       ELSE
 202 *
 203 *        Lower triangular storage: examine D from top to bottom.
 204 *
 205          DO 20 INFO = 1, N
 206             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
 207      $         RETURN
 208    20    CONTINUE
 209       END IF
 210       INFO = 0
 211 *
 212       IF( UPPER ) THEN
 213 *
 214 *        Compute inv(A) from the factorization A = U*D*U**T.
 215 *
 216 *        K is the main loop index, increasing from 1 to N in steps of
 217 *        1 or 2, depending on the size of the diagonal blocks.
 218 *
 219          K = 1
 220    30    CONTINUE
 221 *
 222 *        If K > N, exit from loop.
 223 *
 224          IF( K.GT.N )
 225      $      GO TO 40
 226 *
 227          IF( IPIV( K ).GT.0 ) THEN
 228 *
 229 *           1 x 1 diagonal block
 230 *
 231 *           Invert the diagonal block.
 232 *
 233             A( K, K ) = ONE / A( K, K )
 234 *
 235 *           Compute column K of the inverse.
 236 *
 237             IF( K.GT.1 ) THEN
 238                CALL SCOPY( K-1, A( 1, K ), 1, WORK, 1 )
 239                CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
 240      $                     A( 1, K ), 1 )
 241                A( K, K ) = A( K, K ) - SDOT( K-1, WORK, 1, A( 1, K ),
 242      $                     1 )
 243             END IF
 244             KSTEP = 1
 245          ELSE
 246 *
 247 *           2 x 2 diagonal block
 248 *
 249 *           Invert the diagonal block.
 250 *
 251             T = ABS( A( K, K+1 ) )
 252             AK = A( K, K ) / T
 253             AKP1 = A( K+1, K+1 ) / T
 254             AKKP1 = A( K, K+1 ) / T
 255             D = T*( AK*AKP1-ONE )
 256             A( K, K ) = AKP1 / D
 257             A( K+1, K+1 ) = AK / D
 258             A( K, K+1 ) = -AKKP1 / D
 259 *
 260 *           Compute columns K and K+1 of the inverse.
 261 *
 262             IF( K.GT.1 ) THEN
 263                CALL SCOPY( K-1, A( 1, K ), 1, WORK, 1 )
 264                CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
 265      $                     A( 1, K ), 1 )
 266                A( K, K ) = A( K, K ) - SDOT( K-1, WORK, 1, A( 1, K ),
 267      $                     1 )
 268                A( K, K+1 ) = A( K, K+1 ) -
 269      $                       SDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
 270                CALL SCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
 271                CALL SSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
 272      $                     A( 1, K+1 ), 1 )
 273                A( K+1, K+1 ) = A( K+1, K+1 ) -
 274      $                         SDOT( K-1, WORK, 1, A( 1, K+1 ), 1 )
 275             END IF
 276             KSTEP = 2
 277          END IF
 278 *
 279          IF( KSTEP.EQ.1 ) THEN
 280 *
 281 *           Interchange rows and columns K and IPIV(K) in the leading
 282 *           submatrix A(1:k+1,1:k+1)
 283 *
 284             KP = IPIV( K )
 285             IF( KP.NE.K ) THEN
 286                IF( KP.GT.1 )
 287      $             CALL SSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
 288                CALL SSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
 289                TEMP = A( K, K )
 290                A( K, K ) = A( KP, KP )
 291                A( KP, KP ) = TEMP
 292             END IF
 293          ELSE
 294 *
 295 *           Interchange rows and columns K and K+1 with -IPIV(K) and
 296 *           -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
 297 *
 298             KP = -IPIV( K )
 299             IF( KP.NE.K ) THEN
 300                IF( KP.GT.1 )
 301      $            CALL SSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
 302                CALL SSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
 303 *
 304                TEMP = A( K, K )
 305                A( K, K ) = A( KP, KP )
 306                A( KP, KP ) = TEMP
 307                TEMP = A( K, K+1 )
 308                A( K, K+1 ) = A( KP, K+1 )
 309                A( KP, K+1 ) = TEMP
 310             END IF
 311 *
 312             K = K + 1
 313             KP = -IPIV( K )
 314             IF( KP.NE.K ) THEN
 315                IF( KP.GT.1 )
 316      $            CALL SSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
 317                CALL SSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
 318                TEMP = A( K, K )
 319                A( K, K ) = A( KP, KP )
 320                A( KP, KP ) = TEMP
 321             END IF
 322          END IF
 323 *
 324          K = K + 1
 325          GO TO 30
 326    40    CONTINUE
 327 *
 328       ELSE
 329 *
 330 *        Compute inv(A) from the factorization A = L*D*L**T.
 331 *
 332 *        K is the main loop index, increasing from 1 to N in steps of
 333 *        1 or 2, depending on the size of the diagonal blocks.
 334 *
 335          K = N
 336    50    CONTINUE
 337 *
 338 *        If K < 1, exit from loop.
 339 *
 340          IF( K.LT.1 )
 341      $      GO TO 60
 342 *
 343          IF( IPIV( K ).GT.0 ) THEN
 344 *
 345 *           1 x 1 diagonal block
 346 *
 347 *           Invert the diagonal block.
 348 *
 349             A( K, K ) = ONE / A( K, K )
 350 *
 351 *           Compute column K of the inverse.
 352 *
 353             IF( K.LT.N ) THEN
 354                CALL SCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
 355                CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
 356      $                     ZERO, A( K+1, K ), 1 )
 357                A( K, K ) = A( K, K ) - SDOT( N-K, WORK, 1, A( K+1, K ),
 358      $                     1 )
 359             END IF
 360             KSTEP = 1
 361          ELSE
 362 *
 363 *           2 x 2 diagonal block
 364 *
 365 *           Invert the diagonal block.
 366 *
 367             T = ABS( A( K, K-1 ) )
 368             AK = A( K-1, K-1 ) / T
 369             AKP1 = A( K, K ) / T
 370             AKKP1 = A( K, K-1 ) / T
 371             D = T*( AK*AKP1-ONE )
 372             A( K-1, K-1 ) = AKP1 / D
 373             A( K, K ) = AK / D
 374             A( K, K-1 ) = -AKKP1 / D
 375 *
 376 *           Compute columns K-1 and K of the inverse.
 377 *
 378             IF( K.LT.N ) THEN
 379                CALL SCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
 380                CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
 381      $                     ZERO, A( K+1, K ), 1 )
 382                A( K, K ) = A( K, K ) - SDOT( N-K, WORK, 1, A( K+1, K ),
 383      $                     1 )
 384                A( K, K-1 ) = A( K, K-1 ) -
 385      $                       SDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
 386      $                       1 )
 387                CALL SCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
 388                CALL SSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
 389      $                     ZERO, A( K+1, K-1 ), 1 )
 390                A( K-1, K-1 ) = A( K-1, K-1 ) -
 391      $                         SDOT( N-K, WORK, 1, A( K+1, K-1 ), 1 )
 392             END IF
 393             KSTEP = 2
 394          END IF
 395 *
 396          IF( KSTEP.EQ.1 ) THEN
 397 *
 398 *           Interchange rows and columns K and IPIV(K) in the trailing
 399 *           submatrix A(k-1:n,k-1:n)
 400 *
 401             KP = IPIV( K )
 402             IF( KP.NE.K ) THEN
 403                IF( KP.LT.N )
 404      $            CALL SSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
 405                CALL SSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
 406                TEMP = A( K, K )
 407                A( K, K ) = A( KP, KP )
 408                A( KP, KP ) = TEMP
 409             END IF
 410          ELSE
 411 *
 412 *           Interchange rows and columns K and K-1 with -IPIV(K) and
 413 *           -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
 414 *
 415             KP = -IPIV( K )
 416             IF( KP.NE.K ) THEN
 417                IF( KP.LT.N )
 418      $            CALL SSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
 419                CALL SSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
 420 *
 421                TEMP = A( K, K )
 422                A( K, K ) = A( KP, KP )
 423                A( KP, KP ) = TEMP
 424                TEMP = A( K, K-1 )
 425                A( K, K-1 ) = A( KP, K-1 )
 426                A( KP, K-1 ) = TEMP
 427             END IF
 428 *
 429             K = K - 1
 430             KP = -IPIV( K )
 431             IF( KP.NE.K ) THEN
 432                IF( KP.LT.N )
 433      $            CALL SSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
 434                CALL SSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
 435                TEMP = A( K, K )
 436                A( K, K ) = A( KP, KP )
 437                A( KP, KP ) = TEMP
 438             END IF
 439          END IF
 440 *
 441          K = K - 1
 442          GO TO 50
 443    60    CONTINUE
 444       END IF
 445 *
 446       RETURN
 447 *
 448 *     End of SSYTRI_ROOK
 449 *
 450       END