SRC/dsytf2_rook.f

   1 *> \brief \b DSYTF2_ROOK computes the factorization of a real 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 DSYTF2_ROOK + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2_rook.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2_rook.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2_rook.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSYTF2_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 *       DOUBLE PRECISION   A( LDA, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DSYTF2_ROOK computes the factorization of a real 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 DOUBLE PRECISION 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 doubleSYcomputational
 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 DSYTF2_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       DOUBLE PRECISION   A( LDA, * )
 209 *     ..
 210 *
 211 *  =====================================================================
 212 *
 213 *     .. Parameters ..
 214       DOUBLE PRECISION   ZERO, ONE
 215       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 216       DOUBLE PRECISION   EIGHT, SEVTEN
 217       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
 218 *     ..
 219 *     .. Local Scalars ..
 220       LOGICAL            UPPER, DONE
 221       INTEGER            I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
 222      $                   P, II
 223       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
 224      $                   ROWMAX, DTEMP, T, WK, WKM1, WKP1, SFMIN
 225 *     ..
 226 *     .. External Functions ..
 227       LOGICAL            LSAME
 228       INTEGER            IDAMAX
 229       DOUBLE PRECISION   DLAMCH
 230       EXTERNAL           LSAME, IDAMAX, DLAMCH
 231 *     ..
 232 *     .. External Subroutines ..
 233       EXTERNAL           DSCAL, DSWAP, DSYR, XERBLA
 234 *     ..
 235 *     .. Intrinsic Functions ..
 236       INTRINSIC          ABS, MAX, SQRT
 237 *     ..
 238 *     .. Executable Statements ..
 239 *
 240 *     Test the input parameters.
 241 *
 242       INFO = 0
 243       UPPER = LSAME( UPLO, 'U' )
 244       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 245          INFO = -1
 246       ELSE IF( N.LT.0 ) THEN
 247          INFO = -2
 248       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 249          INFO = -4
 250       END IF
 251       IF( INFO.NE.0 ) THEN
 252          CALL XERBLA( 'DSYTF2_ROOK', -INFO )
 253          RETURN
 254       END IF
 255 *
 256 *     Initialize ALPHA for use in choosing pivot block size.
 257 *
 258       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
 259 *
 260 *     Compute machine safe minimum
 261 *
 262       SFMIN = DLAMCH( 'S' )
 263 *
 264       IF( UPPER ) THEN
 265 *
 266 *        Factorize A as U*D*U**T using the upper triangle of A
 267 *
 268 *        K is the main loop index, decreasing from N to 1 in steps of
 269 *        1 or 2
 270 *
 271          K = N
 272    10    CONTINUE
 273 *
 274 *        If K < 1, exit from loop
 275 *
 276          IF( K.LT.1 )
 277      $      GO TO 70
 278          KSTEP = 1
 279          P = K
 280 *
 281 *        Determine rows and columns to be interchanged and whether
 282 *        a 1-by-1 or 2-by-2 pivot block will be used
 283 *
 284          ABSAKK = ABS( A( K, K ) )
 285 *
 286 *        IMAX is the row-index of the largest off-diagonal element in
 287 *        column K, and COLMAX is its absolute value.
 288 *        Determine both COLMAX and IMAX.
 289 *
 290          IF( K.GT.1 ) THEN
 291             IMAX = IDAMAX( K-1, A( 1, K ), 1 )
 292             COLMAX = ABS( A( IMAX, K ) )
 293          ELSE
 294             COLMAX = ZERO
 295          END IF
 296 *
 297          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN
 298 *
 299 *           Column K is zero or underflow: set INFO and continue
 300 *
 301             IF( INFO.EQ.0 )
 302      $         INFO = K
 303             KP = K
 304          ELSE
 305 *
 306 *           Test for interchange
 307 *
 308 *           Equivalent to testing for (used to handle NaN and Inf)
 309 *           ABSAKK.GE.ALPHA*COLMAX
 310 *
 311             IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
 312 *
 313 *              no interchange,
 314 *              use 1-by-1 pivot block
 315 *
 316                KP = K
 317             ELSE
 318 *
 319                DONE = .FALSE.
 320 *
 321 *              Loop until pivot found
 322 *
 323    12          CONTINUE
 324 *
 325 *                 Begin pivot search loop body
 326 *
 327 *                 JMAX is the column-index of the largest off-diagonal
 328 *                 element in row IMAX, and ROWMAX is its absolute value.
 329 *                 Determine both ROWMAX and JMAX.
 330 *
 331                   IF( IMAX.NE.K ) THEN
 332                      JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ),
 333      $                                    LDA )
 334                      ROWMAX = ABS( A( IMAX, JMAX ) )
 335                   ELSE
 336                      ROWMAX = ZERO
 337                   END IF
 338 *
 339                   IF( IMAX.GT.1 ) THEN
 340                      ITEMP = IDAMAX( IMAX-1, A( 1, IMAX ), 1 )
 341                      DTEMP = ABS( A( ITEMP, IMAX ) )
 342                      IF( DTEMP.GT.ROWMAX ) THEN
 343                         ROWMAX = DTEMP
 344                         JMAX = ITEMP
 345                      END IF
 346                   END IF
 347 *
 348 *                 Equivalent to testing for (used to handle NaN and Inf)
 349 *                 ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
 350 *
 351                   IF( .NOT.( ABS( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ) )
 352      $            THEN
 353 *
 354 *                    interchange rows and columns K and IMAX,
 355 *                    use 1-by-1 pivot block
 356 *
 357                      KP = IMAX
 358                      DONE = .TRUE.
 359 *
 360 *                 Equivalent to testing for ROWMAX .EQ. COLMAX,
 361 *                 used to handle NaN and Inf
 362 *
 363                   ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN
 364 *
 365 *                    interchange rows and columns K+1 and IMAX,
 366 *                    use 2-by-2 pivot block
 367 *
 368                      KP = IMAX
 369                      KSTEP = 2
 370                      DONE = .TRUE.
 371                   ELSE
 372 *
 373 *                    Pivot NOT found, set variables and repeat
 374 *
 375                      P = IMAX
 376                      COLMAX = ROWMAX
 377                      IMAX = JMAX
 378                   END IF
 379 *
 380 *                 End pivot search loop body
 381 *
 382                IF( .NOT. DONE ) GOTO 12
 383 *
 384             END IF
 385 *
 386 *           Swap TWO rows and TWO columns
 387 *
 388 *           First swap
 389 *
 390             IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
 391 *
 392 *              Interchange rows and column K and P in the leading
 393 *              submatrix A(1:k,1:k) if we have a 2-by-2 pivot
 394 *
 395                IF( P.GT.1 )
 396      $            CALL DSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
 397                IF( P.LT.(K-1) )
 398      $            CALL DSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ),
 399      $                     LDA )
 400                T = A( K, K )
 401                A( K, K ) = A( P, P )
 402                A( P, P ) = T
 403             END IF
 404 *
 405 *           Second swap
 406 *
 407             KK = K - KSTEP + 1
 408             IF( KP.NE.KK ) THEN
 409 *
 410 *              Interchange rows and columns KK and KP in the leading
 411 *              submatrix A(1:k,1:k)
 412 *
 413                IF( KP.GT.1 )
 414      $            CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
 415                IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) )
 416      $            CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
 417      $                     LDA )
 418                T = A( KK, KK )
 419                A( KK, KK ) = A( KP, KP )
 420                A( KP, KP ) = T
 421                IF( KSTEP.EQ.2 ) THEN
 422                   T = A( K-1, K )
 423                   A( K-1, K ) = A( KP, K )
 424                   A( KP, K ) = T
 425                END IF
 426             END IF
 427 *
 428 *           Update the leading submatrix
 429 *
 430             IF( KSTEP.EQ.1 ) THEN
 431 *
 432 *              1-by-1 pivot block D(k): column k now holds
 433 *
 434 *              W(k) = U(k)*D(k)
 435 *
 436 *              where U(k) is the k-th column of U
 437 *
 438                IF( K.GT.1 ) THEN
 439 *
 440 *                 Perform a rank-1 update of A(1:k-1,1:k-1) and
 441 *                 store U(k) in column k
 442 *
 443                   IF( ABS( A( K, K ) ).GE.SFMIN ) THEN
 444 *
 445 *                    Perform a rank-1 update of A(1:k-1,1:k-1) as
 446 *                    A := A - U(k)*D(k)*U(k)**T
 447 *                       = A - W(k)*1/D(k)*W(k)**T
 448 *
 449                      D11 = ONE / A( K, K )
 450                      CALL DSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
 451 *
 452 *                    Store U(k) in column k
 453 *
 454                      CALL DSCAL( K-1, D11, A( 1, K ), 1 )
 455                   ELSE
 456 *
 457 *                    Store L(k) in column K
 458 *
 459                      D11 = A( K, K )
 460                      DO 16 II = 1, K - 1
 461                         A( II, K ) = A( II, K ) / D11
 462    16                CONTINUE
 463 *
 464 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 465 *                    A := A - U(k)*D(k)*U(k)**T
 466 *                       = A - W(k)*(1/D(k))*W(k)**T
 467 *                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
 468 *
 469                      CALL DSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
 470                   END IF
 471                END IF
 472 *
 473             ELSE
 474 *
 475 *              2-by-2 pivot block D(k): columns k and k-1 now hold
 476 *
 477 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
 478 *
 479 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
 480 *              of U
 481 *
 482 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
 483 *
 484 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
 485 *                 = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
 486 *
 487 *              and store L(k) and L(k+1) in columns k and k+1
 488 *
 489                IF( K.GT.2 ) THEN
 490 *
 491                   D12 = A( K-1, K )
 492                   D22 = A( K-1, K-1 ) / D12
 493                   D11 = A( K, K ) / D12
 494                   T = ONE / ( D11*D22-ONE )
 495 *
 496                   DO 30 J = K - 2, 1, -1
 497 *
 498                      WKM1 = T*( D11*A( J, K-1 )-A( J, K ) )
 499                      WK = T*( D22*A( J, K )-A( J, K-1 ) )
 500 *
 501                      DO 20 I = J, 1, -1
 502                         A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK -
 503      $                              ( A( I, K-1 ) / D12 )*WKM1
 504    20                CONTINUE
 505 *
 506 *                    Store U(k) and U(k-1) in cols k and k-1 for row J
 507 *
 508                      A( J, K ) = WK / D12
 509                      A( J, K-1 ) = WKM1 / D12
 510 *
 511    30             CONTINUE
 512 *
 513                END IF
 514 *
 515             END IF
 516          END IF
 517 *
 518 *        Store details of the interchanges in IPIV
 519 *
 520          IF( KSTEP.EQ.1 ) THEN
 521             IPIV( K ) = KP
 522          ELSE
 523             IPIV( K ) = -P
 524             IPIV( K-1 ) = -KP
 525          END IF
 526 *
 527 *        Decrease K and return to the start of the main loop
 528 *
 529          K = K - KSTEP
 530          GO TO 10
 531 *
 532       ELSE
 533 *
 534 *        Factorize A as L*D*L**T using the lower triangle of A
 535 *
 536 *        K is the main loop index, increasing from 1 to N in steps of
 537 *        1 or 2
 538 *
 539          K = 1
 540    40    CONTINUE
 541 *
 542 *        If K > N, exit from loop
 543 *
 544          IF( K.GT.N )
 545      $      GO TO 70
 546          KSTEP = 1
 547          P = K
 548 *
 549 *        Determine rows and columns to be interchanged and whether
 550 *        a 1-by-1 or 2-by-2 pivot block will be used
 551 *
 552          ABSAKK = ABS( A( K, K ) )
 553 *
 554 *        IMAX is the row-index of the largest off-diagonal element in
 555 *        column K, and COLMAX is its absolute value.
 556 *        Determine both COLMAX and IMAX.
 557 *
 558          IF( K.LT.N ) THEN
 559             IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 )
 560             COLMAX = ABS( A( IMAX, K ) )
 561          ELSE
 562             COLMAX = ZERO
 563          END IF
 564 *
 565          IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN
 566 *
 567 *           Column K is zero or underflow: set INFO and continue
 568 *
 569             IF( INFO.EQ.0 )
 570      $         INFO = K
 571             KP = K
 572          ELSE
 573 *
 574 *           Test for interchange
 575 *
 576 *           Equivalent to testing for (used to handle NaN and Inf)
 577 *           ABSAKK.GE.ALPHA*COLMAX
 578 *
 579             IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
 580 *
 581 *              no interchange, use 1-by-1 pivot block
 582 *
 583                KP = K
 584             ELSE
 585 *
 586                DONE = .FALSE.
 587 *
 588 *              Loop until pivot found
 589 *
 590    42          CONTINUE
 591 *
 592 *                 Begin pivot search loop body
 593 *
 594 *                 JMAX is the column-index of the largest off-diagonal
 595 *                 element in row IMAX, and ROWMAX is its absolute value.
 596 *                 Determine both ROWMAX and JMAX.
 597 *
 598                   IF( IMAX.NE.K ) THEN
 599                      JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA )
 600                      ROWMAX = ABS( A( IMAX, JMAX ) )
 601                   ELSE
 602                      ROWMAX = ZERO
 603                   END IF
 604 *
 605                   IF( IMAX.LT.N ) THEN
 606                      ITEMP = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ),
 607      $                                     1 )
 608                      DTEMP = ABS( A( ITEMP, IMAX ) )
 609                      IF( DTEMP.GT.ROWMAX ) THEN
 610                         ROWMAX = DTEMP
 611                         JMAX = ITEMP
 612                      END IF
 613                   END IF
 614 *
 615 *                 Equivalent to testing for (used to handle NaN and Inf)
 616 *                 ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
 617 *
 618                   IF( .NOT.( ABS( A( IMAX, IMAX ) ).LT.ALPHA*ROWMAX ) )
 619      $            THEN
 620 *
 621 *                    interchange rows and columns K and IMAX,
 622 *                    use 1-by-1 pivot block
 623 *
 624                      KP = IMAX
 625                      DONE = .TRUE.
 626 *
 627 *                 Equivalent to testing for ROWMAX .EQ. COLMAX,
 628 *                 used to handle NaN and Inf
 629 *
 630                   ELSE IF( ( P.EQ.JMAX ).OR.( ROWMAX.LE.COLMAX ) ) THEN
 631 *
 632 *                    interchange rows and columns K+1 and IMAX,
 633 *                    use 2-by-2 pivot block
 634 *
 635                      KP = IMAX
 636                      KSTEP = 2
 637                      DONE = .TRUE.
 638                   ELSE
 639 *
 640 *                    Pivot NOT found, set variables and repeat
 641 *
 642                      P = IMAX
 643                      COLMAX = ROWMAX
 644                      IMAX = JMAX
 645                   END IF
 646 *
 647 *                 End pivot search loop body
 648 *
 649                IF( .NOT. DONE ) GOTO 42
 650 *
 651             END IF
 652 *
 653 *           Swap TWO rows and TWO columns
 654 *
 655 *           First swap
 656 *
 657             IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
 658 *
 659 *              Interchange rows and column K and P in the trailing
 660 *              submatrix A(k:n,k:n) if we have a 2-by-2 pivot
 661 *
 662                IF( P.LT.N )
 663      $            CALL DSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
 664                IF( P.GT.(K+1) )
 665      $            CALL DSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA )
 666                T = A( K, K )
 667                A( K, K ) = A( P, P )
 668                A( P, P ) = T
 669             END IF
 670 *
 671 *           Second swap
 672 *
 673             KK = K + KSTEP - 1
 674             IF( KP.NE.KK ) THEN
 675 *
 676 *              Interchange rows and columns KK and KP in the trailing
 677 *              submatrix A(k:n,k:n)
 678 *
 679                IF( KP.LT.N )
 680      $            CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
 681                IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) )
 682      $            CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
 683      $                     LDA )
 684                T = A( KK, KK )
 685                A( KK, KK ) = A( KP, KP )
 686                A( KP, KP ) = T
 687                IF( KSTEP.EQ.2 ) THEN
 688                   T = A( K+1, K )
 689                   A( K+1, K ) = A( KP, K )
 690                   A( KP, K ) = T
 691                END IF
 692             END IF
 693 *
 694 *           Update the trailing submatrix
 695 *
 696             IF( KSTEP.EQ.1 ) THEN
 697 *
 698 *              1-by-1 pivot block D(k): column k now holds
 699 *
 700 *              W(k) = L(k)*D(k)
 701 *
 702 *              where L(k) is the k-th column of L
 703 *
 704                IF( K.LT.N ) THEN
 705 *
 706 *              Perform a rank-1 update of A(k+1:n,k+1:n) and
 707 *              store L(k) in column k
 708 *
 709                   IF( ABS( A( K, K ) ).GE.SFMIN ) THEN
 710 *
 711 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 712 *                    A := A - L(k)*D(k)*L(k)**T
 713 *                       = A - W(k)*(1/D(k))*W(k)**T
 714 *
 715                      D11 = ONE / A( K, K )
 716                      CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
 717      $                          A( K+1, K+1 ), LDA )
 718 *
 719 *                    Store L(k) in column k
 720 *
 721                      CALL DSCAL( N-K, D11, A( K+1, K ), 1 )
 722                   ELSE
 723 *
 724 *                    Store L(k) in column k
 725 *
 726                      D11 = A( K, K )
 727                      DO 46 II = K + 1, N
 728                         A( II, K ) = A( II, K ) / D11
 729    46                CONTINUE
 730 *
 731 *                    Perform a rank-1 update of A(k+1:n,k+1:n) as
 732 *                    A := A - L(k)*D(k)*L(k)**T
 733 *                       = A - W(k)*(1/D(k))*W(k)**T
 734 *                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
 735 *
 736                      CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
 737      $                          A( K+1, K+1 ), LDA )
 738                   END IF
 739                END IF
 740 *
 741             ELSE
 742 *
 743 *              2-by-2 pivot block D(k): columns k and k+1 now hold
 744 *
 745 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
 746 *
 747 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
 748 *              of L
 749 *
 750 *
 751 *              Perform a rank-2 update of A(k+2:n,k+2:n) as
 752 *
 753 *              A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
 754 *                 = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
 755 *
 756 *              and store L(k) and L(k+1) in columns k and k+1
 757 *
 758                IF( K.LT.N-1 ) THEN
 759 *
 760                   D21 = A( K+1, K )
 761                   D11 = A( K+1, K+1 ) / D21
 762                   D22 = A( K, K ) / D21
 763                   T = ONE / ( D11*D22-ONE )
 764 *
 765                   DO 60 J = K + 2, N
 766 *
 767 *                    Compute  D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
 768 *
 769                      WK = T*( D11*A( J, K )-A( J, K+1 ) )
 770                      WKP1 = T*( D22*A( J, K+1 )-A( J, K ) )
 771 *
 772 *                    Perform a rank-2 update of A(k+2:n,k+2:n)
 773 *
 774                      DO 50 I = J, N
 775                         A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK -
 776      $                              ( A( I, K+1 ) / D21 )*WKP1
 777    50                CONTINUE
 778 *
 779 *                    Store L(k) and L(k+1) in cols k and k+1 for row J
 780 *
 781                      A( J, K ) = WK / D21
 782                      A( J, K+1 ) = WKP1 / D21
 783 *
 784    60             CONTINUE
 785 *
 786                END IF
 787 *
 788             END IF
 789          END IF
 790 *
 791 *        Store details of the interchanges in IPIV
 792 *
 793          IF( KSTEP.EQ.1 ) THEN
 794             IPIV( K ) = KP
 795          ELSE
 796             IPIV( K ) = -P
 797             IPIV( K+1 ) = -KP
 798          END IF
 799 *
 800 *        Increase K and return to the start of the main loop
 801 *
 802          K = K + KSTEP
 803          GO TO 40
 804 *
 805       END IF
 806 *
 807    70 CONTINUE
 808 *
 809       RETURN
 810 *
 811 *     End of DSYTF2_ROOK
 812 *
 813       END