SRC/chetf2_rook.f

   1 *> \brief \b CHETF2_ROOK computes the factorization of a complex Hermitian 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 CHETF2_ROOK + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetf2_rook.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetf2_rook.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetf2_rook.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHETF2_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 *> CHETF2_ROOK computes the factorization of a complex Hermitian matrix A
  39 *> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:
  40 *>
  41 *>    A = U*D*U**H  or  A = L*D*L**H
  42 *>
  43 *> where U (or L) is a product of permutation and unit upper (lower)
  44 *> triangular matrices, U**H is the conjugate transpose of U, and D is
  45 *> Hermitian and 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 *>          Hermitian 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 Hermitian 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) < 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 complexHEcomputational
 135 *
 136 *> \par Further Details:
 137 *  =====================
 138 *>
 139 *> \verbatim
 140 *>
 141 *>  If UPLO = 'U', then A = U*D*U**H, 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**H, 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, USA
 192 *> \endverbatim
 193 *
 194 *  =====================================================================
 195       SUBROUTINE CHETF2_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 *     ..
 219 *     .. Local Scalars ..
 220       LOGICAL            DONE, UPPER
 221       INTEGER            I, II, IMAX, ITEMP, J, JMAX, K, KK, KP, KSTEP,
 222      $                   P
 223       REAL               ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, STEMP,
 224      $                   ROWMAX, TT, SFMIN
 225       COMPLEX            D12, D21, T, WK, WKM1, WKP1, Z
 226 *     ..
 227 *     .. External Functions ..
 228 *
 229       LOGICAL            LSAME
 230       INTEGER            ICAMAX
 231       REAL               SLAMCH, SLAPY2
 232       EXTERNAL           LSAME, ICAMAX, SLAMCH, SLAPY2
 233 *     ..
 234 *     .. External Subroutines ..
 235       EXTERNAL           XERBLA, CSSCAL, CHER, CSWAP
 236 *     ..
 237 *     .. Intrinsic Functions ..
 238       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL, SQRT
 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( 'CHETF2_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**H 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 = ABS( REAL( 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             A( K, K ) = REAL( A( K, K ) )
 313          ELSE
 314 *
 315 *           ============================================================
 316 *
 317 *           BEGIN pivot search
 318 *
 319 *           Case(1)
 320 *           Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
 321 *           (used to handle NaN and Inf)
 322 *
 323             IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
 324 *
 325 *              no interchange, use 1-by-1 pivot block
 326 *
 327                KP = K
 328 *
 329             ELSE
 330 *
 331                DONE = .FALSE.
 332 *
 333 *              Loop until pivot found
 334 *
 335    12          CONTINUE
 336 *
 337 *                 BEGIN pivot search loop body
 338 *
 339 *
 340 *                 JMAX is the column-index of the largest off-diagonal
 341 *                 element in row IMAX, and ROWMAX is its absolute value.
 342 *                 Determine both ROWMAX and JMAX.
 343 *
 344                   IF( IMAX.NE.K ) THEN
 345                      JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ),
 346      $                                     LDA )
 347                      ROWMAX = CABS1( A( IMAX, JMAX ) )
 348                   ELSE
 349                      ROWMAX = ZERO
 350                   END IF
 351 *
 352                   IF( IMAX.GT.1 ) THEN
 353                      ITEMP = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
 354                      STEMP = CABS1( A( ITEMP, IMAX ) )
 355                      IF( STEMP.GT.ROWMAX ) THEN
 356                         ROWMAX = STEMP
 357                         JMAX = ITEMP
 358                      END IF
 359                   END IF
 360 *
 361 *                 Case(2)
 362 *                 Equivalent to testing for
 363 *                 ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
 364 *                 (used to handle NaN and Inf)
 365 *
 366                   IF( .NOT.( ABS( REAL( A( IMAX, IMAX ) ) )
 367      $                       .LT.ALPHA*ROWMAX ) ) THEN
 368 *
 369 *                    interchange rows and columns K and IMAX,
 370 *                    use 1-by-1 pivot block
 371 *
 372                      KP = IMAX
 373                      DONE = .TRUE.
 374 *
 375 *                 Case(3)
 376 *                 Equivalent to testing for ROWMAX.EQ.COLMAX,
 377 *                 (used to handle NaN and Inf)
 378 *
 379                   ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
 380      $            THEN
 381 *
 382 *                    interchange rows and columns K-1 and IMAX,
 383 *                    use 2-by-2 pivot block
 384 *
 385                      KP = IMAX
 386                      KSTEP = 2
 387                      DONE = .TRUE.
 388 *
 389 *                 Case(4)
 390                   ELSE
 391 *
 392 *                    Pivot not found: set params and repeat
 393 *
 394                      P = IMAX
 395                      COLMAX = ROWMAX
 396                      IMAX = JMAX
 397                   END IF
 398 *
 399 *                 END pivot search loop body
 400 *
 401                IF( .NOT.DONE ) GOTO 12
 402 *
 403             END IF
 404 *
 405 *           END pivot search
 406 *
 407 *           ============================================================
 408 *
 409 *           KK is the column of A where pivoting step stopped
 410 *
 411             KK = K - KSTEP + 1
 412 *
 413 *           For only a 2x2 pivot, interchange rows and columns K and P
 414 *           in the leading submatrix A(1:k,1:k)
 415 *
 416             IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
 417 *              (1) Swap columnar parts
 418                IF( P.GT.1 )
 419      $            CALL CSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
 420 *              (2) Swap and conjugate middle parts
 421                DO 14 J = P + 1, K - 1
 422                   T = CONJG( A( J, K ) )
 423                   A( J, K ) = CONJG( A( P, J ) )
 424                   A( P, J ) = T
 425    14          CONTINUE
 426 *              (3) Swap and conjugate corner elements at row-col interserction
 427                A( P, K ) = CONJG( A( P, K ) )
 428 *              (4) Swap diagonal elements at row-col intersection
 429                R1 = REAL( A( K, K ) )
 430                A( K, K ) = REAL( A( P, P ) )
 431                A( P, P ) = R1
 432             END IF
 433 *
 434 *           For both 1x1 and 2x2 pivots, interchange rows and
 435 *           columns KK and KP in the leading submatrix A(1:k,1:k)
 436 *
 437             IF( KP.NE.KK ) THEN
 438 *              (1) Swap columnar parts
 439                IF( KP.GT.1 )
 440      $            CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
 441 *              (2) Swap and conjugate middle parts
 442                DO 15 J = KP + 1, KK - 1
 443                   T = CONJG( A( J, KK ) )
 444                   A( J, KK ) = CONJG( A( KP, J ) )
 445                   A( KP, J ) = T
 446    15          CONTINUE
 447 *              (3) Swap and conjugate corner elements at row-col interserction
 448                A( KP, KK ) = CONJG( A( KP, KK ) )
 449 *              (4) Swap diagonal elements at row-col intersection
 450                R1 = REAL( A( KK, KK ) )
 451                A( KK, KK ) = REAL( A( KP, KP ) )
 452                A( KP, KP ) = R1
 453 *
 454                IF( KSTEP.EQ.2 ) THEN
 455 *                 (*) Make sure that diagonal element of pivot is real
 456                   A( K, K ) = REAL( A( K, K ) )
 457 *                 (5) Swap row elements
 458                   T = A( K-1, K )
 459                   A( K-1, K ) = A( KP, K )
 460                   A( KP, K ) = T
 461                END IF
 462             ELSE
 463 *              (*) Make sure that diagonal element of pivot is real
 464                A( K, K ) = REAL( A( K, K ) )
 465                IF( KSTEP.EQ.2 )
 466      $            A( K-1, K-1 ) = REAL( A( K-1, K-1 ) )
 467             END IF
 468 *
 469 *           Update the leading submatrix
 470 *
 471             IF( KSTEP.EQ.1 ) THEN
 472 *
 473 *              1-by-1 pivot block D(k): column k now holds
 474 *
 475 *              W(k) = U(k)*D(k)
 476 *
 477 *              where U(k) is the k-th column of U
 478 *
 479                IF( K.GT.1 ) THEN
 480 *
 481 *                 Perform a rank-1 update of A(1:k-1,1:k-1) and
 482 *                 store U(k) in column k
 483 *
 484                   IF( ABS( REAL( A( K, K ) ) ).GE.SFMIN ) THEN
 485 *
 486 *                    Perform a rank-1 update of A(1:k-1,1:k-1) as
 487 *                    A := A - U(k)*D(k)*U(k)**T
 488 *                       = A - W(k)*1/D(k)*W(k)**T
 489 *
 490                      D11 = ONE / REAL( A( K, K ) )
 491                      CALL CHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
 492 *
 493 *                    Store U(k) in column k
 494 *
 495                      CALL CSSCAL( K-1, D11, A( 1, K ), 1 )
 496                   ELSE
 497 *
 498 *                    Store L(k) in column K
 499 *
 500                      D11 = REAL( A( K, K ) )
 501                      DO 16 II = 1, K - 1
 502                         A( II, K ) = A( II, K ) / D11
 503    16                CONTINUE
 504 *
 505 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 506 *                    A := A - U(k)*D(k)*U(k)**T
 507 *                       = A - W(k)*(1/D(k))*W(k)**T
 508 *                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
 509 *
 510                      CALL CHER( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
 511                   END IF
 512                END IF
 513 *
 514             ELSE
 515 *
 516 *              2-by-2 pivot block D(k): columns k and k-1 now hold
 517 *
 518 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
 519 *
 520 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
 521 *              of U
 522 *
 523 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
 524 *
 525 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
 526 *                 = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
 527 *
 528 *              and store L(k) and L(k+1) in columns k and k+1
 529 *
 530                IF( K.GT.2 ) THEN
 531 *                 D = |A12|
 532                   D = SLAPY2( REAL( A( K-1, K ) ),
 533      $                AIMAG( A( K-1, K ) ) )
 534                   D11 = A( K, K ) / D
 535                   D22 = A( K-1, K-1 ) / D
 536                   D12 = A( K-1, K ) / D
 537                   TT = ONE / ( D11*D22-ONE )
 538 *
 539                   DO 30 J = K - 2, 1, -1
 540 *
 541 *                    Compute  D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
 542 *
 543                      WKM1 = TT*( D11*A( J, K-1 )-CONJG( D12 )*
 544      $                      A( J, K ) )
 545                      WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) )
 546 *
 547 *                    Perform a rank-2 update of A(1:k-2,1:k-2)
 548 *
 549                      DO 20 I = J, 1, -1
 550                         A( I, J ) = A( I, J ) -
 551      $                              ( A( I, K ) / D )*CONJG( WK ) -
 552      $                              ( A( I, K-1 ) / D )*CONJG( WKM1 )
 553    20                CONTINUE
 554 *
 555 *                    Store U(k) and U(k-1) in cols k and k-1 for row J
 556 *
 557                      A( J, K ) = WK / D
 558                      A( J, K-1 ) = WKM1 / D
 559 *                    (*) Make sure that diagonal element of pivot is real
 560                      A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO )
 561 *
 562    30             CONTINUE
 563 *
 564                END IF
 565 *
 566             END IF
 567 *
 568          END IF
 569 *
 570 *        Store details of the interchanges in IPIV
 571 *
 572          IF( KSTEP.EQ.1 ) THEN
 573             IPIV( K ) = KP
 574          ELSE
 575             IPIV( K ) = -P
 576             IPIV( K-1 ) = -KP
 577          END IF
 578 *
 579 *        Decrease K and return to the start of the main loop
 580 *
 581          K = K - KSTEP
 582          GO TO 10
 583 *
 584       ELSE
 585 *
 586 *        Factorize A as L*D*L**H using the lower triangle of A
 587 *
 588 *        K is the main loop index, increasing from 1 to N in steps of
 589 *        1 or 2
 590 *
 591          K = 1
 592    40    CONTINUE
 593 *
 594 *        If K > N, exit from loop
 595 *
 596          IF( K.GT.N )
 597      $      GO TO 70
 598          KSTEP = 1
 599          P = K
 600 *
 601 *        Determine rows and columns to be interchanged and whether
 602 *        a 1-by-1 or 2-by-2 pivot block will be used
 603 *
 604          ABSAKK = ABS( REAL( A( K, K ) ) )
 605 *
 606 *        IMAX is the row-index of the largest off-diagonal element in
 607 *        column K, and COLMAX is its absolute value.
 608 *        Determine both COLMAX and IMAX.
 609 *
 610          IF( K.LT.N ) THEN
 611             IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
 612             COLMAX = CABS1( A( IMAX, K ) )
 613          ELSE
 614             COLMAX = ZERO
 615          END IF
 616 *
 617          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
 618 *
 619 *           Column K is zero or underflow: set INFO and continue
 620 *
 621             IF( INFO.EQ.0 )
 622      $         INFO = K
 623             KP = K
 624             A( K, K ) = REAL( A( K, K ) )
 625          ELSE
 626 *
 627 *           ============================================================
 628 *
 629 *           BEGIN pivot search
 630 *
 631 *           Case(1)
 632 *           Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
 633 *           (used to handle NaN and Inf)
 634 *
 635             IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
 636 *
 637 *              no interchange, use 1-by-1 pivot block
 638 *
 639                KP = K
 640 *
 641             ELSE
 642 *
 643                DONE = .FALSE.
 644 *
 645 *              Loop until pivot found
 646 *
 647    42          CONTINUE
 648 *
 649 *                 BEGIN pivot search loop body
 650 *
 651 *
 652 *                 JMAX is the column-index of the largest off-diagonal
 653 *                 element in row IMAX, and ROWMAX is its absolute value.
 654 *                 Determine both ROWMAX and JMAX.
 655 *
 656                   IF( IMAX.NE.K ) THEN
 657                      JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
 658                      ROWMAX = CABS1( A( IMAX, JMAX ) )
 659                   ELSE
 660                      ROWMAX = ZERO
 661                   END IF
 662 *
 663                   IF( IMAX.LT.N ) THEN
 664                      ITEMP = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ),
 665      $                                     1 )
 666                      STEMP = CABS1( A( ITEMP, IMAX ) )
 667                      IF( STEMP.GT.ROWMAX ) THEN
 668                         ROWMAX = STEMP
 669                         JMAX = ITEMP
 670                      END IF
 671                   END IF
 672 *
 673 *                 Case(2)
 674 *                 Equivalent to testing for
 675 *                 ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
 676 *                 (used to handle NaN and Inf)
 677 *
 678                   IF( .NOT.( ABS( REAL( A( IMAX, IMAX ) ) )
 679      $                       .LT.ALPHA*ROWMAX ) ) THEN
 680 *
 681 *                    interchange rows and columns K and IMAX,
 682 *                    use 1-by-1 pivot block
 683 *
 684                      KP = IMAX
 685                      DONE = .TRUE.
 686 *
 687 *                 Case(3)
 688 *                 Equivalent to testing for ROWMAX.EQ.COLMAX,
 689 *                 (used to handle NaN and Inf)
 690 *
 691                   ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
 692      $            THEN
 693 *
 694 *                    interchange rows and columns K+1 and IMAX,
 695 *                    use 2-by-2 pivot block
 696 *
 697                      KP = IMAX
 698                      KSTEP = 2
 699                      DONE = .TRUE.
 700 *
 701 *                 Case(4)
 702                   ELSE
 703 *
 704 *                    Pivot not found: set params and repeat
 705 *
 706                      P = IMAX
 707                      COLMAX = ROWMAX
 708                      IMAX = JMAX
 709                   END IF
 710 *
 711 *
 712 *                 END pivot search loop body
 713 *
 714                IF( .NOT.DONE ) GOTO 42
 715 *
 716             END IF
 717 *
 718 *           END pivot search
 719 *
 720 *           ============================================================
 721 *
 722 *           KK is the column of A where pivoting step stopped
 723 *
 724             KK = K + KSTEP - 1
 725 *
 726 *           For only a 2x2 pivot, interchange rows and columns K and P
 727 *           in the trailing submatrix A(k:n,k:n)
 728 *
 729             IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
 730 *              (1) Swap columnar parts
 731                IF( P.LT.N )
 732      $            CALL CSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
 733 *              (2) Swap and conjugate middle parts
 734                DO 44 J = K + 1, P - 1
 735                   T = CONJG( A( J, K ) )
 736                   A( J, K ) = CONJG( A( P, J ) )
 737                   A( P, J ) = T
 738    44          CONTINUE
 739 *              (3) Swap and conjugate corner elements at row-col interserction
 740                A( P, K ) = CONJG( A( P, K ) )
 741 *              (4) Swap diagonal elements at row-col intersection
 742                R1 = REAL( A( K, K ) )
 743                A( K, K ) = REAL( A( P, P ) )
 744                A( P, P ) = R1
 745             END IF
 746 *
 747 *           For both 1x1 and 2x2 pivots, interchange rows and
 748 *           columns KK and KP in the trailing submatrix A(k:n,k:n)
 749 *
 750             IF( KP.NE.KK ) THEN
 751 *              (1) Swap columnar parts
 752                IF( KP.LT.N )
 753      $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
 754 *              (2) Swap and conjugate middle parts
 755                DO 45 J = KK + 1, KP - 1
 756                   T = CONJG( A( J, KK ) )
 757                   A( J, KK ) = CONJG( A( KP, J ) )
 758                   A( KP, J ) = T
 759    45          CONTINUE
 760 *              (3) Swap and conjugate corner elements at row-col interserction
 761                A( KP, KK ) = CONJG( A( KP, KK ) )
 762 *              (4) Swap diagonal elements at row-col intersection
 763                R1 = REAL( A( KK, KK ) )
 764                A( KK, KK ) = REAL( A( KP, KP ) )
 765                A( KP, KP ) = R1
 766 *
 767                IF( KSTEP.EQ.2 ) THEN
 768 *                 (*) Make sure that diagonal element of pivot is real
 769                   A( K, K ) = REAL( A( K, K ) )
 770 *                 (5) Swap row elements
 771                   T = A( K+1, K )
 772                   A( K+1, K ) = A( KP, K )
 773                   A( KP, K ) = T
 774                END IF
 775             ELSE
 776 *              (*) Make sure that diagonal element of pivot is real
 777                A( K, K ) = REAL( A( K, K ) )
 778                IF( KSTEP.EQ.2 )
 779      $            A( K+1, K+1 ) = REAL( A( K+1, K+1 ) )
 780             END IF
 781 *
 782 *           Update the trailing submatrix
 783 *
 784             IF( KSTEP.EQ.1 ) THEN
 785 *
 786 *              1-by-1 pivot block D(k): column k of A now holds
 787 *
 788 *              W(k) = L(k)*D(k),
 789 *
 790 *              where L(k) is the k-th column of L
 791 *
 792                IF( K.LT.N ) THEN
 793 *
 794 *                 Perform a rank-1 update of A(k+1:n,k+1:n) and
 795 *                 store L(k) in column k
 796 *
 797 *                 Handle division by a small number
 798 *
 799                   IF( ABS( REAL( A( K, K ) ) ).GE.SFMIN ) THEN
 800 *
 801 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 802 *                    A := A - L(k)*D(k)*L(k)**T
 803 *                       = A - W(k)*(1/D(k))*W(k)**T
 804 *
 805                      D11 = ONE / REAL( A( K, K ) )
 806                      CALL CHER( UPLO, N-K, -D11, A( K+1, K ), 1,
 807      $                          A( K+1, K+1 ), LDA )
 808 *
 809 *                    Store L(k) in column k
 810 *
 811                      CALL CSSCAL( N-K, D11, A( K+1, K ), 1 )
 812                   ELSE
 813 *
 814 *                    Store L(k) in column k
 815 *
 816                      D11 = REAL( A( K, K ) )
 817                      DO 46 II = K + 1, N
 818                         A( II, K ) = A( II, K ) / D11
 819    46                CONTINUE
 820 *
 821 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 822 *                    A := A - L(k)*D(k)*L(k)**T
 823 *                       = A - W(k)*(1/D(k))*W(k)**T
 824 *                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
 825 *
 826                      CALL CHER( UPLO, N-K, -D11, A( K+1, K ), 1,
 827      $                          A( K+1, K+1 ), LDA )
 828                   END IF
 829                END IF
 830 *
 831             ELSE
 832 *
 833 *              2-by-2 pivot block D(k): columns k and k+1 now hold
 834 *
 835 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
 836 *
 837 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
 838 *              of L
 839 *
 840 *
 841 *              Perform a rank-2 update of A(k+2:n,k+2:n) as
 842 *
 843 *              A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
 844 *                 = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
 845 *
 846 *              and store L(k) and L(k+1) in columns k and k+1
 847 *
 848                IF( K.LT.N-1 ) THEN
 849 *                 D = |A21|
 850                   D = SLAPY2( REAL( A( K+1, K ) ),
 851      $                AIMAG( A( K+1, K ) ) )
 852                   D11 = REAL( A( K+1, K+1 ) ) / D
 853                   D22 = REAL( A( K, K ) ) / D
 854                   D21 = A( K+1, K ) / D
 855                   TT = ONE / ( D11*D22-ONE )
 856 *
 857                   DO 60 J = K + 2, N
 858 *
 859 *                    Compute  D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
 860 *
 861                      WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) )
 862                      WKP1 = TT*( D22*A( J, K+1 )-CONJG( D21 )*
 863      $                      A( J, K ) )
 864 *
 865 *                    Perform a rank-2 update of A(k+2:n,k+2:n)
 866 *
 867                      DO 50 I = J, N
 868                         A( I, J ) = A( I, J ) -
 869      $                              ( A( I, K ) / D )*CONJG( WK ) -
 870      $                              ( A( I, K+1 ) / D )*CONJG( WKP1 )
 871    50                CONTINUE
 872 *
 873 *                    Store L(k) and L(k+1) in cols k and k+1 for row J
 874 *
 875                      A( J, K ) = WK / D
 876                      A( J, K+1 ) = WKP1 / D
 877 *                    (*) Make sure that diagonal element of pivot is real
 878                      A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO )
 879 *
 880    60             CONTINUE
 881 *
 882                END IF
 883 *
 884             END IF
 885 *
 886          END IF
 887 *
 888 *        Store details of the interchanges in IPIV
 889 *
 890          IF( KSTEP.EQ.1 ) THEN
 891             IPIV( K ) = KP
 892          ELSE
 893             IPIV( K ) = -P
 894             IPIV( K+1 ) = -KP
 895          END IF
 896 *
 897 *        Increase K and return to the start of the main loop
 898 *
 899          K = K + KSTEP
 900          GO TO 40
 901 *
 902       END IF
 903 *
 904    70 CONTINUE
 905 *
 906       RETURN
 907 *
 908 *     End of CHETF2_ROOK
 909 *
 910       END