SRC/csytf2_rook.f

   1 *> \brief \b CSYTF2_ROOK computes the factorization of a complex symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") 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_ROOK + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2_rook.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2_rook.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2_rook.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CSYTF2_ROOK( 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_ROOK computes the factorization of a complex symmetric matrix A
  39 *> using the bounded Bunch-Kaufman ("rook") 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)
  96 *>             were interchanged and D(k,k) is a 1-by-1 diagonal block.
  97 *>
  98 *>             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
  99 *>             columns k and -IPIV(k) were interchanged and rows and
 100 *>             columns k-1 and -IPIV(k-1) were inerchaged,
 101 *>             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
 102 *>
 103 *>          If UPLO = 'L':
 104 *>             If IPIV(k) > 0, then rows and columns k and IPIV(k)
 105 *>             were interchanged and D(k,k) is a 1-by-1 diagonal block.
 106 *>
 107 *>             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
 108 *>             columns k and -IPIV(k) were interchanged and rows and
 109 *>             columns k+1 and -IPIV(k+1) were inerchaged,
 110 *>             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
 111 *> \endverbatim
 112 *>
 113 *> \param[out] INFO
 114 *> \verbatim
 115 *>          INFO is INTEGER
 116 *>          = 0: successful exit
 117 *>          < 0: if INFO = -k, the k-th argument had an illegal value
 118 *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
 119 *>               has been completed, but the block diagonal matrix D is
 120 *>               exactly singular, and division by zero will occur if it
 121 *>               is used to solve a system of equations.
 122 *> \endverbatim
 123 *
 124 *  Authors:
 125 *  ========
 126 *
 127 *> \author Univ. of Tennessee
 128 *> \author Univ. of California Berkeley
 129 *> \author Univ. of Colorado Denver
 130 *> \author NAG Ltd.
 131 *
 132 *> \date November 2013
 133 *
 134 *> \ingroup complexSYcomputational
 135 *
 136 *> \par Further Details:
 137 *  =====================
 138 *>
 139 *> \verbatim
 140 *>
 141 *>  If UPLO = 'U', then A = U*D*U**T, where
 142 *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
 143 *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
 144 *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
 145 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
 146 *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
 147 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
 148 *>
 149 *>             (   I    v    0   )   k-s
 150 *>     U(k) =  (   0    I    0   )   s
 151 *>             (   0    0    I   )   n-k
 152 *>                k-s   s   n-k
 153 *>
 154 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
 155 *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
 156 *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
 157 *>
 158 *>  If UPLO = 'L', then A = L*D*L**T, where
 159 *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
 160 *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
 161 *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
 162 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
 163 *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
 164 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
 165 *>
 166 *>             (   I    0     0   )  k-1
 167 *>     L(k) =  (   0    I     0   )  s
 168 *>             (   0    v     I   )  n-k-s+1
 169 *>                k-1   s  n-k-s+1
 170 *>
 171 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
 172 *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
 173 *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
 174 *> \endverbatim
 175 *
 176 *> \par Contributors:
 177 *  ==================
 178 *>
 179 *> \verbatim
 180 *>
 181 *>  November 2013,     Igor Kozachenko,
 182 *>                  Computer Science Division,
 183 *>                  University of California, Berkeley
 184 *>
 185 *>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
 186 *>                  School of Mathematics,
 187 *>                  University of Manchester
 188 *>
 189 *>  01-01-96 - Based on modifications by
 190 *>    J. Lewis, Boeing Computer Services Company
 191 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA
 192 *> \endverbatim
 193 *
 194 *  =====================================================================
 195       SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO )
 196 *
 197 *  -- LAPACK computational routine (version 3.5.0) --
 198 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 199 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 200 *     November 2013
 201 *
 202 *     .. Scalar Arguments ..
 203       CHARACTER          UPLO
 204       INTEGER            INFO, LDA, N
 205 *     ..
 206 *     .. Array Arguments ..
 207       INTEGER            IPIV( * )
 208       COMPLEX            A( LDA, * )
 209 *     ..
 210 *
 211 *  =====================================================================
 212 *
 213 *     .. Parameters ..
 214       REAL               ZERO, ONE
 215       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 216       REAL               EIGHT, SEVTEN
 217       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
 218       COMPLEX            CONE
 219       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 220 *     ..
 221 *     .. Local Scalars ..
 222       LOGICAL            UPPER, DONE
 223       INTEGER            I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
 224      $                   P, II
 225       REAL               ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
 226       COMPLEX            D11, D12, D21, D22, T, WK, WKM1, WKP1, Z
 227 *     ..
 228 *     .. External Functions ..
 229       LOGICAL            LSAME
 230       INTEGER            ICAMAX
 231       REAL               SLAMCH
 232       EXTERNAL           LSAME, ICAMAX, SLAMCH
 233 *     ..
 234 *     .. External Subroutines ..
 235       EXTERNAL           CSCAL, CSWAP, CSYR, XERBLA
 236 *     ..
 237 *     .. Intrinsic Functions ..
 238       INTRINSIC          ABS, MAX, SQRT, AIMAG, REAL
 239 *     ..
 240 *     .. Statement Functions ..
 241       REAL               CABS1
 242 *     ..
 243 *     .. Statement Function definitions ..
 244       CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
 245 *     ..
 246 *     .. Executable Statements ..
 247 *
 248 *     Test the input parameters.
 249 *
 250       INFO = 0
 251       UPPER = LSAME( UPLO, 'U' )
 252       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 253          INFO = -1
 254       ELSE IF( N.LT.0 ) THEN
 255          INFO = -2
 256       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 257          INFO = -4
 258       END IF
 259       IF( INFO.NE.0 ) THEN
 260          CALL XERBLA( 'CSYTF2_ROOK', -INFO )
 261          RETURN
 262       END IF
 263 *
 264 *     Initialize ALPHA for use in choosing pivot block size.
 265 *
 266       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
 267 *
 268 *     Compute machine safe minimum
 269 *
 270       SFMIN = SLAMCH( 'S' )
 271 *
 272       IF( UPPER ) THEN
 273 *
 274 *        Factorize A as U*D*U**T using the upper triangle of A
 275 *
 276 *        K is the main loop index, decreasing from N to 1 in steps of
 277 *        1 or 2
 278 *
 279          K = N
 280    10    CONTINUE
 281 *
 282 *        If K < 1, exit from loop
 283 *
 284          IF( K.LT.1 )
 285      $      GO TO 70
 286          KSTEP = 1
 287          P = K
 288 *
 289 *        Determine rows and columns to be interchanged and whether
 290 *        a 1-by-1 or 2-by-2 pivot block will be used
 291 *
 292          ABSAKK = CABS1( A( K, K ) )
 293 *
 294 *        IMAX is the row-index of the largest off-diagonal element in
 295 *        column K, and COLMAX is its absolute value.
 296 *        Determine both COLMAX and IMAX.
 297 *
 298          IF( K.GT.1 ) THEN
 299             IMAX = ICAMAX( K-1, A( 1, K ), 1 )
 300             COLMAX = CABS1( A( IMAX, K ) )
 301          ELSE
 302             COLMAX = ZERO
 303          END IF
 304 *
 305          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN
 306 *
 307 *           Column K is zero or underflow: set INFO and continue
 308 *
 309             IF( INFO.EQ.0 )
 310      $         INFO = K
 311             KP = K
 312          ELSE
 313 *
 314 *           Test for interchange
 315 *
 316 *           Equivalent to testing for (used to handle NaN and Inf)
 317 *           ABSAKK.GE.ALPHA*COLMAX
 318 *
 319             IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
 320 *
 321 *              no interchange,
 322 *              use 1-by-1 pivot block
 323 *
 324                KP = K
 325             ELSE
 326 *
 327                DONE = .FALSE.
 328 *
 329 *              Loop until pivot found
 330 *
 331    12          CONTINUE
 332 *
 333 *                 Begin pivot search loop body
 334 *
 335 *                 JMAX is the column-index of the largest off-diagonal
 336 *                 element in row IMAX, and ROWMAX is its absolute value.
 337 *                 Determine both ROWMAX and JMAX.
 338 *
 339                   IF( IMAX.NE.K ) THEN
 340                      JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ),
 341      $                                    LDA )
 342                      ROWMAX = CABS1( A( IMAX, JMAX ) )
 343                   ELSE
 344                      ROWMAX = ZERO
 345                   END IF
 346 *
 347                   IF( IMAX.GT.1 ) THEN
 348                      ITEMP = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
 349                      STEMP = CABS1( A( ITEMP, IMAX ) )
 350                      IF( STEMP.GT.ROWMAX ) THEN
 351                         ROWMAX = STEMP
 352                         JMAX = ITEMP
 353                      END IF
 354                   END IF
 355 *
 356 *                 Equivalent to testing for (used to handle NaN and Inf)
 357 *                 CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
 358 *
 359                   IF( .NOT.( CABS1(A( IMAX, IMAX )).LT.ALPHA*ROWMAX ) )
 360      $            THEN
 361 *
 362 *                    interchange rows and columns K and IMAX,
 363 *                    use 1-by-1 pivot block
 364 *
 365                      KP = IMAX
 366                      DONE = .TRUE.
 367 *
 368 *                 Equivalent to testing for ROWMAX .EQ. COLMAX,
 369 *                 used to handle NaN and Inf
 370 *
 371                   ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN
 372 *
 373 *                    interchange rows and columns K+1 and IMAX,
 374 *                    use 2-by-2 pivot block
 375 *
 376                      KP = IMAX
 377                      KSTEP = 2
 378                      DONE = .TRUE.
 379                   ELSE
 380 *
 381 *                    Pivot NOT found, set variables and repeat
 382 *
 383                      P = IMAX
 384                      COLMAX = ROWMAX
 385                      IMAX = JMAX
 386                   END IF
 387 *
 388 *                 End pivot search loop body
 389 *
 390                IF( .NOT. DONE ) GOTO 12
 391 *
 392             END IF
 393 *
 394 *           Swap TWO rows and TWO columns
 395 *
 396 *           First swap
 397 *
 398             IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
 399 *
 400 *              Interchange rows and column K and P in the leading
 401 *              submatrix A(1:k,1:k) if we have a 2-by-2 pivot
 402 *
 403                IF( P.GT.1 )
 404      $            CALL CSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
 405                IF( P.LT.(K-1) )
 406      $            CALL CSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ),
 407      $                     LDA )
 408                T = A( K, K )
 409                A( K, K ) = A( P, P )
 410                A( P, P ) = T
 411             END IF
 412 *
 413 *           Second swap
 414 *
 415             KK = K - KSTEP + 1
 416             IF( KP.NE.KK ) THEN
 417 *
 418 *              Interchange rows and columns KK and KP in the leading
 419 *              submatrix A(1:k,1:k)
 420 *
 421                IF( KP.GT.1 )
 422      $            CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
 423                IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) )
 424      $            CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
 425      $                     LDA )
 426                T = A( KK, KK )
 427                A( KK, KK ) = A( KP, KP )
 428                A( KP, KP ) = T
 429                IF( KSTEP.EQ.2 ) THEN
 430                   T = A( K-1, K )
 431                   A( K-1, K ) = A( KP, K )
 432                   A( KP, K ) = T
 433                END IF
 434             END IF
 435 *
 436 *           Update the leading submatrix
 437 *
 438             IF( KSTEP.EQ.1 ) THEN
 439 *
 440 *              1-by-1 pivot block D(k): column k now holds
 441 *
 442 *              W(k) = U(k)*D(k)
 443 *
 444 *              where U(k) is the k-th column of U
 445 *
 446                IF( K.GT.1 ) THEN
 447 *
 448 *                 Perform a rank-1 update of A(1:k-1,1:k-1) and
 449 *                 store U(k) in column k
 450 *
 451                   IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
 452 *
 453 *                    Perform a rank-1 update of A(1:k-1,1:k-1) as
 454 *                    A := A - U(k)*D(k)*U(k)**T
 455 *                       = A - W(k)*1/D(k)*W(k)**T
 456 *
 457                      D11 = CONE / A( K, K )
 458                      CALL CSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
 459 *
 460 *                    Store U(k) in column k
 461 *
 462                      CALL CSCAL( K-1, D11, A( 1, K ), 1 )
 463                   ELSE
 464 *
 465 *                    Store L(k) in column K
 466 *
 467                      D11 = A( K, K )
 468                      DO 16 II = 1, K - 1
 469                         A( II, K ) = A( II, K ) / D11
 470    16                CONTINUE
 471 *
 472 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 473 *                    A := A - U(k)*D(k)*U(k)**T
 474 *                       = A - W(k)*(1/D(k))*W(k)**T
 475 *                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
 476 *
 477                      CALL CSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
 478                   END IF
 479                END IF
 480 *
 481             ELSE
 482 *
 483 *              2-by-2 pivot block D(k): columns k and k-1 now hold
 484 *
 485 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
 486 *
 487 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
 488 *              of U
 489 *
 490 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
 491 *
 492 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
 493 *                 = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
 494 *
 495 *              and store L(k) and L(k+1) in columns k and k+1
 496 *
 497                IF( K.GT.2 ) THEN
 498 *
 499                   D12 = A( K-1, K )
 500                   D22 = A( K-1, K-1 ) / D12
 501                   D11 = A( K, K ) / D12
 502                   T = CONE / ( D11*D22-CONE )
 503 *
 504                   DO 30 J = K - 2, 1, -1
 505 *
 506                      WKM1 = T*( D11*A( J, K-1 )-A( J, K ) )
 507                      WK = T*( D22*A( J, K )-A( J, K-1 ) )
 508 *
 509                      DO 20 I = J, 1, -1
 510                         A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK -
 511      $                              ( A( I, K-1 ) / D12 )*WKM1
 512    20                CONTINUE
 513 *
 514 *                    Store U(k) and U(k-1) in cols k and k-1 for row J
 515 *
 516                      A( J, K ) = WK / D12
 517                      A( J, K-1 ) = WKM1 / D12
 518 *
 519    30             CONTINUE
 520 *
 521                END IF
 522 *
 523             END IF
 524          END IF
 525 *
 526 *        Store details of the interchanges in IPIV
 527 *
 528          IF( KSTEP.EQ.1 ) THEN
 529             IPIV( K ) = KP
 530          ELSE
 531             IPIV( K ) = -P
 532             IPIV( K-1 ) = -KP
 533          END IF
 534 *
 535 *        Decrease K and return to the start of the main loop
 536 *
 537          K = K - KSTEP
 538          GO TO 10
 539 *
 540       ELSE
 541 *
 542 *        Factorize A as L*D*L**T using the lower triangle of A
 543 *
 544 *        K is the main loop index, increasing from 1 to N in steps of
 545 *        1 or 2
 546 *
 547          K = 1
 548    40    CONTINUE
 549 *
 550 *        If K > N, exit from loop
 551 *
 552          IF( K.GT.N )
 553      $      GO TO 70
 554          KSTEP = 1
 555          P = K
 556 *
 557 *        Determine rows and columns to be interchanged and whether
 558 *        a 1-by-1 or 2-by-2 pivot block will be used
 559 *
 560          ABSAKK = CABS1( A( K, K ) )
 561 *
 562 *        IMAX is the row-index of the largest off-diagonal element in
 563 *        column K, and COLMAX is its absolute value.
 564 *        Determine both COLMAX and IMAX.
 565 *
 566          IF( K.LT.N ) THEN
 567             IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
 568             COLMAX = CABS1( A( IMAX, K ) )
 569          ELSE
 570             COLMAX = ZERO
 571          END IF
 572 *
 573          IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN
 574 *
 575 *           Column K is zero or underflow: set INFO and continue
 576 *
 577             IF( INFO.EQ.0 )
 578      $         INFO = K
 579             KP = K
 580          ELSE
 581 *
 582 *           Test for interchange
 583 *
 584 *           Equivalent to testing for (used to handle NaN and Inf)
 585 *           ABSAKK.GE.ALPHA*COLMAX
 586 *
 587             IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
 588 *
 589 *              no interchange, use 1-by-1 pivot block
 590 *
 591                KP = K
 592             ELSE
 593 *
 594                DONE = .FALSE.
 595 *
 596 *              Loop until pivot found
 597 *
 598    42          CONTINUE
 599 *
 600 *                 Begin pivot search loop body
 601 *
 602 *                 JMAX is the column-index of the largest off-diagonal
 603 *                 element in row IMAX, and ROWMAX is its absolute value.
 604 *                 Determine both ROWMAX and JMAX.
 605 *
 606                   IF( IMAX.NE.K ) THEN
 607                      JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
 608                      ROWMAX = CABS1( A( IMAX, JMAX ) )
 609                   ELSE
 610                      ROWMAX = ZERO
 611                   END IF
 612 *
 613                   IF( IMAX.LT.N ) THEN
 614                      ITEMP = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ),
 615      $                                     1 )
 616                      STEMP = CABS1( A( ITEMP, IMAX ) )
 617                      IF( STEMP.GT.ROWMAX ) THEN
 618                         ROWMAX = STEMP
 619                         JMAX = ITEMP
 620                      END IF
 621                   END IF
 622 *
 623 *                 Equivalent to testing for (used to handle NaN and Inf)
 624 *                 CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
 625 *
 626                   IF( .NOT.( CABS1(A( IMAX, IMAX )).LT.ALPHA*ROWMAX ) )
 627      $            THEN
 628 *
 629 *                    interchange rows and columns K and IMAX,
 630 *                    use 1-by-1 pivot block
 631 *
 632                      KP = IMAX
 633                      DONE = .TRUE.
 634 *
 635 *                 Equivalent to testing for ROWMAX .EQ. COLMAX,
 636 *                 used to handle NaN and Inf
 637 *
 638                   ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN
 639 *
 640 *                    interchange rows and columns K+1 and IMAX,
 641 *                    use 2-by-2 pivot block
 642 *
 643                      KP = IMAX
 644                      KSTEP = 2
 645                      DONE = .TRUE.
 646                   ELSE
 647 *
 648 *                    Pivot NOT found, set variables and repeat
 649 *
 650                      P = IMAX
 651                      COLMAX = ROWMAX
 652                      IMAX = JMAX
 653                   END IF
 654 *
 655 *                 End pivot search loop body
 656 *
 657                IF( .NOT. DONE ) GOTO 42
 658 *
 659             END IF
 660 *
 661 *           Swap TWO rows and TWO columns
 662 *
 663 *           First swap
 664 *
 665             IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
 666 *
 667 *              Interchange rows and column K and P in the trailing
 668 *              submatrix A(k:n,k:n) if we have a 2-by-2 pivot
 669 *
 670                IF( P.LT.N )
 671      $            CALL CSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
 672                IF( P.GT.(K+1) )
 673      $            CALL CSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA )
 674                T = A( K, K )
 675                A( K, K ) = A( P, P )
 676                A( P, P ) = T
 677             END IF
 678 *
 679 *           Second swap
 680 *
 681             KK = K + KSTEP - 1
 682             IF( KP.NE.KK ) THEN
 683 *
 684 *              Interchange rows and columns KK and KP in the trailing
 685 *              submatrix A(k:n,k:n)
 686 *
 687                IF( KP.LT.N )
 688      $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
 689                IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) )
 690      $            CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
 691      $                     LDA )
 692                T = A( KK, KK )
 693                A( KK, KK ) = A( KP, KP )
 694                A( KP, KP ) = T
 695                IF( KSTEP.EQ.2 ) THEN
 696                   T = A( K+1, K )
 697                   A( K+1, K ) = A( KP, K )
 698                   A( KP, K ) = T
 699                END IF
 700             END IF
 701 *
 702 *           Update the trailing submatrix
 703 *
 704             IF( KSTEP.EQ.1 ) THEN
 705 *
 706 *              1-by-1 pivot block D(k): column k now holds
 707 *
 708 *              W(k) = L(k)*D(k)
 709 *
 710 *              where L(k) is the k-th column of L
 711 *
 712                IF( K.LT.N ) THEN
 713 *
 714 *              Perform a rank-1 update of A(k+1:n,k+1:n) and
 715 *              store L(k) in column k
 716 *
 717                   IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
 718 *
 719 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 720 *                    A := A - L(k)*D(k)*L(k)**T
 721 *                       = A - W(k)*(1/D(k))*W(k)**T
 722 *
 723                      D11 = CONE / A( K, K )
 724                      CALL CSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
 725      $                          A( K+1, K+1 ), LDA )
 726 *
 727 *                    Store L(k) in column k
 728 *
 729                      CALL CSCAL( N-K, D11, A( K+1, K ), 1 )
 730                   ELSE
 731 *
 732 *                    Store L(k) in column k
 733 *
 734                      D11 = A( K, K )
 735                      DO 46 II = K + 1, N
 736                         A( II, K ) = A( II, K ) / D11
 737    46                CONTINUE
 738 *
 739 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 740 *                    A := A - L(k)*D(k)*L(k)**T
 741 *                       = A - W(k)*(1/D(k))*W(k)**T
 742 *                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
 743 *
 744                      CALL CSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
 745      $                          A( K+1, K+1 ), LDA )
 746                   END IF
 747                END IF
 748 *
 749             ELSE
 750 *
 751 *              2-by-2 pivot block D(k): columns k and k+1 now hold
 752 *
 753 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
 754 *
 755 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
 756 *              of L
 757 *
 758 *
 759 *              Perform a rank-2 update of A(k+2:n,k+2:n) as
 760 *
 761 *              A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
 762 *                 = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
 763 *
 764 *              and store L(k) and L(k+1) in columns k and k+1
 765 *
 766                IF( K.LT.N-1 ) THEN
 767 *
 768                   D21 = A( K+1, K )
 769                   D11 = A( K+1, K+1 ) / D21
 770                   D22 = A( K, K ) / D21
 771                   T = CONE / ( D11*D22-CONE )
 772 *
 773                   DO 60 J = K + 2, N
 774 *
 775 *                    Compute  D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
 776 *
 777                      WK = T*( D11*A( J, K )-A( J, K+1 ) )
 778                      WKP1 = T*( D22*A( J, K+1 )-A( J, K ) )
 779 *
 780 *                    Perform a rank-2 update of A(k+2:n,k+2:n)
 781 *
 782                      DO 50 I = J, N
 783                         A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK -
 784      $                              ( A( I, K+1 ) / D21 )*WKP1
 785    50                CONTINUE
 786 *
 787 *                    Store L(k) and L(k+1) in cols k and k+1 for row J
 788 *
 789                      A( J, K ) = WK / D21
 790                      A( J, K+1 ) = WKP1 / D21
 791 *
 792    60             CONTINUE
 793 *
 794                END IF
 795 *
 796             END IF
 797          END IF
 798 *
 799 *        Store details of the interchanges in IPIV
 800 *
 801          IF( KSTEP.EQ.1 ) THEN
 802             IPIV( K ) = KP
 803          ELSE
 804             IPIV( K ) = -P
 805             IPIV( K+1 ) = -KP
 806          END IF
 807 *
 808 *        Increase K and return to the start of the main loop
 809 *
 810          K = K + KSTEP
 811          GO TO 40
 812 *
 813       END IF
 814 *
 815    70 CONTINUE
 816 *
 817       RETURN
 818 *
 819 *     End of CSYTF2_ROOK
 820 *
 821       END