SRC/zsytri_rook.f

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