SRC/csytf2.f

   1 *> \brief \b CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CSYTF2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, LDA, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IPIV( * )
  29 *       COMPLEX            A( LDA, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CSYTF2 computes the factorization of a complex symmetric matrix A
  39 *> using the Bunch-Kaufman diagonal pivoting method:
  40 *>
  41 *>    A = U*D*U**T  or  A = L*D*L**T
  42 *>
  43 *> where U (or L) is a product of permutation and unit upper (lower)
  44 *> triangular matrices, U**T is the transpose of U, and D is symmetric and
  45 *> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
  46 *>
  47 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
  48 *> \endverbatim
  49 *
  50 *  Arguments:
  51 *  ==========
  52 *
  53 *> \param[in] UPLO
  54 *> \verbatim
  55 *>          UPLO is CHARACTER*1
  56 *>          Specifies whether the upper or lower triangular part of the
  57 *>          symmetric matrix A is stored:
  58 *>          = 'U':  Upper triangular
  59 *>          = 'L':  Lower triangular
  60 *> \endverbatim
  61 *>
  62 *> \param[in] N
  63 *> \verbatim
  64 *>          N is INTEGER
  65 *>          The order of the matrix A.  N >= 0.
  66 *> \endverbatim
  67 *>
  68 *> \param[in,out] A
  69 *> \verbatim
  70 *>          A is COMPLEX array, dimension (LDA,N)
  71 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  72 *>          n-by-n upper triangular part of A contains the upper
  73 *>          triangular part of the matrix A, and the strictly lower
  74 *>          triangular part of A is not referenced.  If UPLO = 'L', the
  75 *>          leading n-by-n lower triangular part of A contains the lower
  76 *>          triangular part of the matrix A, and the strictly upper
  77 *>          triangular part of A is not referenced.
  78 *>
  79 *>          On exit, the block diagonal matrix D and the multipliers used
  80 *>          to obtain the factor U or L (see below for further details).
  81 *> \endverbatim
  82 *>
  83 *> \param[in] LDA
  84 *> \verbatim
  85 *>          LDA is INTEGER
  86 *>          The leading dimension of the array A.  LDA >= max(1,N).
  87 *> \endverbatim
  88 *>
  89 *> \param[out] IPIV
  90 *> \verbatim
  91 *>          IPIV is INTEGER array, dimension (N)
  92 *>          Details of the interchanges and the block structure of D.
  93 *>
  94 *>          If UPLO = 'U':
  95 *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  96 *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
  97 *>
  98 *>             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
  99 *>             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
 100 *>             is a 2-by-2 diagonal block.
 101 *>
 102 *>          If UPLO = 'L':
 103 *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
 104 *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
 105 *>
 106 *>             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
 107 *>             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
 108 *>             is a 2-by-2 diagonal block.
 109 *> \endverbatim
 110 *>
 111 *> \param[out] INFO
 112 *> \verbatim
 113 *>          INFO is INTEGER
 114 *>          = 0: successful exit
 115 *>          < 0: if INFO = -k, the k-th argument had an illegal value
 116 *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
 117 *>               has been completed, but the block diagonal matrix D is
 118 *>               exactly singular, and division by zero will occur if it
 119 *>               is used to solve a system of equations.
 120 *> \endverbatim
 121 *
 122 *  Authors:
 123 *  ========
 124 *
 125 *> \author Univ. of Tennessee
 126 *> \author Univ. of California Berkeley
 127 *> \author Univ. of Colorado Denver
 128 *> \author NAG Ltd.
 129 *
 130 *> \date November 2013
 131 *
 132 *> \ingroup complexSYcomputational
 133 *
 134 *> \par Further Details:
 135 *  =====================
 136 *>
 137 *> \verbatim
 138 *>
 139 *>  If UPLO = 'U', then A = U*D*U**T, where
 140 *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
 141 *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
 142 *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
 143 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
 144 *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
 145 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
 146 *>
 147 *>             (   I    v    0   )   k-s
 148 *>     U(k) =  (   0    I    0   )   s
 149 *>             (   0    0    I   )   n-k
 150 *>                k-s   s   n-k
 151 *>
 152 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
 153 *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
 154 *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
 155 *>
 156 *>  If UPLO = 'L', then A = L*D*L**T, where
 157 *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
 158 *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
 159 *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
 160 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
 161 *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
 162 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
 163 *>
 164 *>             (   I    0     0   )  k-1
 165 *>     L(k) =  (   0    I     0   )  s
 166 *>             (   0    v     I   )  n-k-s+1
 167 *>                k-1   s  n-k-s+1
 168 *>
 169 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
 170 *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
 171 *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
 172 *> \endverbatim
 173 *
 174 *> \par Contributors:
 175 *  ==================
 176 *>
 177 *> \verbatim
 178 *>
 179 *>  09-29-06 - patch from
 180 *>    Bobby Cheng, MathWorks
 181 *>
 182 *>    Replace l.209 and l.377
 183 *>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
 184 *>    by
 185 *>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
 186 *>
 187 *>  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
 188 *>         Company
 189 *> \endverbatim
 190 *
 191 *  =====================================================================
 192       SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO )
 193 *
 194 *  -- LAPACK computational routine (version 3.5.0) --
 195 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 196 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 197 *     November 2013
 198 *
 199 *     .. Scalar Arguments ..
 200       CHARACTER          UPLO
 201       INTEGER            INFO, LDA, N
 202 *     ..
 203 *     .. Array Arguments ..
 204       INTEGER            IPIV( * )
 205       COMPLEX            A( LDA, * )
 206 *     ..
 207 *
 208 *  =====================================================================
 209 *
 210 *     .. Parameters ..
 211       REAL               ZERO, ONE
 212       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 213       REAL               EIGHT, SEVTEN
 214       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
 215       COMPLEX            CONE
 216       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 217 *     ..
 218 *     .. Local Scalars ..
 219       LOGICAL            UPPER
 220       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
 221       REAL               ABSAKK, ALPHA, COLMAX, ROWMAX
 222       COMPLEX            D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
 223 *     ..
 224 *     .. External Functions ..
 225       LOGICAL            LSAME, SISNAN
 226       INTEGER            ICAMAX
 227       EXTERNAL           LSAME, ICAMAX, SISNAN
 228 *     ..
 229 *     .. External Subroutines ..
 230       EXTERNAL           CSCAL, CSWAP, CSYR, XERBLA
 231 *     ..
 232 *     .. Intrinsic Functions ..
 233       INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
 234 *     ..
 235 *     .. Statement Functions ..
 236       REAL               CABS1
 237 *     ..
 238 *     .. Statement Function definitions ..
 239       CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
 240 *     ..
 241 *     .. Executable Statements ..
 242 *
 243 *     Test the input parameters.
 244 *
 245       INFO = 0
 246       UPPER = LSAME( UPLO, 'U' )
 247       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 248          INFO = -1
 249       ELSE IF( N.LT.0 ) THEN
 250          INFO = -2
 251       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 252          INFO = -4
 253       END IF
 254       IF( INFO.NE.0 ) THEN
 255          CALL XERBLA( 'CSYTF2', -INFO )
 256          RETURN
 257       END IF
 258 *
 259 *     Initialize ALPHA for use in choosing pivot block size.
 260 *
 261       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
 262 *
 263       IF( UPPER ) THEN
 264 *
 265 *        Factorize A as U*D*U**T using the upper triangle of A
 266 *
 267 *        K is the main loop index, decreasing from N to 1 in steps of
 268 *        1 or 2
 269 *
 270          K = N
 271    10    CONTINUE
 272 *
 273 *        If K < 1, exit from loop
 274 *
 275          IF( K.LT.1 )
 276      $      GO TO 70
 277          KSTEP = 1
 278 *
 279 *        Determine rows and columns to be interchanged and whether
 280 *        a 1-by-1 or 2-by-2 pivot block will be used
 281 *
 282          ABSAKK = CABS1( A( K, K ) )
 283 *
 284 *        IMAX is the row-index of the largest off-diagonal element in
 285 *        column K, and COLMAX is its absolute value.
 286 *        Determine both COLMAX and IMAX.
 287 *
 288          IF( K.GT.1 ) THEN
 289             IMAX = ICAMAX( K-1, A( 1, K ), 1 )
 290             COLMAX = CABS1( A( IMAX, K ) )
 291          ELSE
 292             COLMAX = ZERO
 293          END IF
 294 *
 295          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN
 296 *
 297 *           Column K is zero or underflow, or contains a NaN:
 298 *           set INFO and continue
 299 *
 300             IF( INFO.EQ.0 )
 301      $         INFO = K
 302             KP = K
 303          ELSE
 304             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
 305 *
 306 *              no interchange, use 1-by-1 pivot block
 307 *
 308                KP = K
 309             ELSE
 310 *
 311 *              JMAX is the column-index of the largest off-diagonal
 312 *              element in row IMAX, and ROWMAX is its absolute value
 313 *
 314                JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
 315                ROWMAX = CABS1( A( IMAX, JMAX ) )
 316                IF( IMAX.GT.1 ) THEN
 317                   JMAX = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
 318                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
 319                END IF
 320 *
 321                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
 322 *
 323 *                 no interchange, use 1-by-1 pivot block
 324 *
 325                   KP = K
 326                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
 327 *
 328 *                 interchange rows and columns K and IMAX, use 1-by-1
 329 *                 pivot block
 330 *
 331                   KP = IMAX
 332                ELSE
 333 *
 334 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
 335 *                 pivot block
 336 *
 337                   KP = IMAX
 338                   KSTEP = 2
 339                END IF
 340             END IF
 341 *
 342             KK = K - KSTEP + 1
 343             IF( KP.NE.KK ) THEN
 344 *
 345 *              Interchange rows and columns KK and KP in the leading
 346 *              submatrix A(1:k,1:k)
 347 *
 348                CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
 349                CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
 350      $                     LDA )
 351                T = A( KK, KK )
 352                A( KK, KK ) = A( KP, KP )
 353                A( KP, KP ) = T
 354                IF( KSTEP.EQ.2 ) THEN
 355                   T = A( K-1, K )
 356                   A( K-1, K ) = A( KP, K )
 357                   A( KP, K ) = T
 358                END IF
 359             END IF
 360 *
 361 *           Update the leading submatrix
 362 *
 363             IF( KSTEP.EQ.1 ) THEN
 364 *
 365 *              1-by-1 pivot block D(k): column k now holds
 366 *
 367 *              W(k) = U(k)*D(k)
 368 *
 369 *              where U(k) is the k-th column of U
 370 *
 371 *              Perform a rank-1 update of A(1:k-1,1:k-1) as
 372 *
 373 *              A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
 374 *
 375                R1 = CONE / A( K, K )
 376                CALL CSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
 377 *
 378 *              Store U(k) in column k
 379 *
 380                CALL CSCAL( K-1, R1, A( 1, K ), 1 )
 381             ELSE
 382 *
 383 *              2-by-2 pivot block D(k): columns k and k-1 now hold
 384 *
 385 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
 386 *
 387 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
 388 *              of U
 389 *
 390 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
 391 *
 392 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
 393 *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
 394 *
 395                IF( K.GT.2 ) THEN
 396 *
 397                   D12 = A( K-1, K )
 398                   D22 = A( K-1, K-1 ) / D12
 399                   D11 = A( K, K ) / D12
 400                   T = CONE / ( D11*D22-CONE )
 401                   D12 = T / D12
 402 *
 403                   DO 30 J = K - 2, 1, -1
 404                      WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
 405                      WK = D12*( D22*A( J, K )-A( J, K-1 ) )
 406                      DO 20 I = J, 1, -1
 407                         A( I, J ) = A( I, J ) - A( I, K )*WK -
 408      $                              A( I, K-1 )*WKM1
 409    20                CONTINUE
 410                      A( J, K ) = WK
 411                      A( J, K-1 ) = WKM1
 412    30             CONTINUE
 413 *
 414                END IF
 415 *
 416             END IF
 417          END IF
 418 *
 419 *        Store details of the interchanges in IPIV
 420 *
 421          IF( KSTEP.EQ.1 ) THEN
 422             IPIV( K ) = KP
 423          ELSE
 424             IPIV( K ) = -KP
 425             IPIV( K-1 ) = -KP
 426          END IF
 427 *
 428 *        Decrease K and return to the start of the main loop
 429 *
 430          K = K - KSTEP
 431          GO TO 10
 432 *
 433       ELSE
 434 *
 435 *        Factorize A as L*D*L**T using the lower triangle of A
 436 *
 437 *        K is the main loop index, increasing from 1 to N in steps of
 438 *        1 or 2
 439 *
 440          K = 1
 441    40    CONTINUE
 442 *
 443 *        If K > N, exit from loop
 444 *
 445          IF( K.GT.N )
 446      $      GO TO 70
 447          KSTEP = 1
 448 *
 449 *        Determine rows and columns to be interchanged and whether
 450 *        a 1-by-1 or 2-by-2 pivot block will be used
 451 *
 452          ABSAKK = CABS1( A( K, K ) )
 453 *
 454 *        IMAX is the row-index of the largest off-diagonal element in
 455 *        column K, and COLMAX is its absolute value.
 456 *        Determine both COLMAX and IMAX.
 457 *
 458          IF( K.LT.N ) THEN
 459             IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
 460             COLMAX = CABS1( A( IMAX, K ) )
 461          ELSE
 462             COLMAX = ZERO
 463          END IF
 464 *
 465          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN
 466 *
 467 *           Column K is zero or underflow, or contains a NaN:
 468 *           set INFO and continue
 469 *
 470             IF( INFO.EQ.0 )
 471      $         INFO = K
 472             KP = K
 473          ELSE
 474             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
 475 *
 476 *              no interchange, use 1-by-1 pivot block
 477 *
 478                KP = K
 479             ELSE
 480 *
 481 *              JMAX is the column-index of the largest off-diagonal
 482 *              element in row IMAX, and ROWMAX is its absolute value
 483 *
 484                JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
 485                ROWMAX = CABS1( A( IMAX, JMAX ) )
 486                IF( IMAX.LT.N ) THEN
 487                   JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
 488                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
 489                END IF
 490 *
 491                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
 492 *
 493 *                 no interchange, use 1-by-1 pivot block
 494 *
 495                   KP = K
 496                ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
 497 *
 498 *                 interchange rows and columns K and IMAX, use 1-by-1
 499 *                 pivot block
 500 *
 501                   KP = IMAX
 502                ELSE
 503 *
 504 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
 505 *                 pivot block
 506 *
 507                   KP = IMAX
 508                   KSTEP = 2
 509                END IF
 510             END IF
 511 *
 512             KK = K + KSTEP - 1
 513             IF( KP.NE.KK ) THEN
 514 *
 515 *              Interchange rows and columns KK and KP in the trailing
 516 *              submatrix A(k:n,k:n)
 517 *
 518                IF( KP.LT.N )
 519      $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
 520                CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
 521      $                     LDA )
 522                T = A( KK, KK )
 523                A( KK, KK ) = A( KP, KP )
 524                A( KP, KP ) = T
 525                IF( KSTEP.EQ.2 ) THEN
 526                   T = A( K+1, K )
 527                   A( K+1, K ) = A( KP, K )
 528                   A( KP, K ) = T
 529                END IF
 530             END IF
 531 *
 532 *           Update the trailing submatrix
 533 *
 534             IF( KSTEP.EQ.1 ) THEN
 535 *
 536 *              1-by-1 pivot block D(k): column k now holds
 537 *
 538 *              W(k) = L(k)*D(k)
 539 *
 540 *              where L(k) is the k-th column of L
 541 *
 542                IF( K.LT.N ) THEN
 543 *
 544 *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
 545 *
 546 *                 A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
 547 *
 548                   R1 = CONE / A( K, K )
 549                   CALL CSYR( UPLO, N-K, -R1, A( K+1, K ), 1,
 550      $                       A( K+1, K+1 ), LDA )
 551 *
 552 *                 Store L(k) in column K
 553 *
 554                   CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
 555                END IF
 556             ELSE
 557 *
 558 *              2-by-2 pivot block D(k)
 559 *
 560                IF( K.LT.N-1 ) THEN
 561 *
 562 *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
 563 *
 564 *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
 565 *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
 566 *
 567 *                 where L(k) and L(k+1) are the k-th and (k+1)-th
 568 *                 columns of L
 569 *
 570                   D21 = A( K+1, K )
 571                   D11 = A( K+1, K+1 ) / D21
 572                   D22 = A( K, K ) / D21
 573                   T = CONE / ( D11*D22-CONE )
 574                   D21 = T / D21
 575 *
 576                   DO 60 J = K + 2, N
 577                      WK = D21*( D11*A( J, K )-A( J, K+1 ) )
 578                      WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
 579                      DO 50 I = J, N
 580                         A( I, J ) = A( I, J ) - A( I, K )*WK -
 581      $                              A( I, K+1 )*WKP1
 582    50                CONTINUE
 583                      A( J, K ) = WK
 584                      A( J, K+1 ) = WKP1
 585    60             CONTINUE
 586                END IF
 587             END IF
 588          END IF
 589 *
 590 *        Store details of the interchanges in IPIV
 591 *
 592          IF( KSTEP.EQ.1 ) THEN
 593             IPIV( K ) = KP
 594          ELSE
 595             IPIV( K ) = -KP
 596             IPIV( K+1 ) = -KP
 597          END IF
 598 *
 599 *        Increase K and return to the start of the main loop
 600 *
 601          K = K + KSTEP
 602          GO TO 40
 603 *
 604       END IF
 605 *
 606    70 CONTINUE
 607       RETURN
 608 *
 609 *     End of CSYTF2
 610 *
 611       END