text: *> \brief \b CSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm).
  *
  *  =========== DOCUMENTATION ===========
  *
  * Online html documentation available at
  *            http://www.netlib.org/lapack/explore-html/
  *
  *> \htmlonly
  *> Download CSYTF2 + dependencies
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2.f">
  *> [TGZ]</a>
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2.f">
  *> [ZIP]</a>
  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2.f">
  *> [TXT]</a>
  *> \endhtmlonly
  *
  *  Definition:
  *  ===========
  *
  *       SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO )
  *
  *       .. Scalar Arguments ..
  *       CHARACTER          UPLO
  *       INTEGER            INFO, LDA, N
  *       ..
  *       .. Array Arguments ..
  *       INTEGER            IPIV( * )
  *       COMPLEX            A( LDA, * )
  *       ..
  *
  *
  *> \par Purpose:
  *  =============
  *>
  *> \verbatim
  *>
  *> CSYTF2 computes the factorization of a complex symmetric matrix A
  *> using the Bunch-Kaufman diagonal pivoting method:
  *>
  *>    A = U*D*U**T  or  A = L*D*L**T
  *>
  *> where U (or L) is a product of permutation and unit upper (lower)
  *> triangular matrices, U**T is the transpose of U, and D is symmetric and
  *> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
  *>
  *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
  *> \endverbatim
  *
  *  Arguments:
  *  ==========
  *
  *> \param[in] UPLO
  *> \verbatim
  *>          UPLO is CHARACTER*1
  *>          Specifies whether the upper or lower triangular part of the
  *>          symmetric matrix A is stored:
  *>          = 'U':  Upper triangular
  *>          = 'L':  Lower triangular
  *> \endverbatim
  *>
  *> \param[in] N
  *> \verbatim
  *>          N is INTEGER
  *>          The order of the matrix A.  N >= 0.
  *> \endverbatim
  *>
  *> \param[in,out] A
  *> \verbatim
  *>          A is COMPLEX array, dimension (LDA,N)
  *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  *>          n-by-n upper triangular part of A contains the upper
  *>          triangular part of the matrix A, and the strictly lower
  *>          triangular part of A is not referenced.  If UPLO = 'L', the
  *>          leading n-by-n lower triangular part of A contains the lower
  *>          triangular part of the matrix A, and the strictly upper
  *>          triangular part of A is not referenced.
  *>
  *>          On exit, the block diagonal matrix D and the multipliers used
  *>          to obtain the factor U or L (see below for further details).
  *> \endverbatim
  *>
  *> \param[in] LDA
  *> \verbatim
  *>          LDA is INTEGER
  *>          The leading dimension of the array A.  LDA >= max(1,N).
  *> \endverbatim
  *>
  *> \param[out] IPIV
  *> \verbatim
  *>          IPIV is INTEGER array, dimension (N)
  *>          Details of the interchanges and the block structure of D.
  *>
  *>          If UPLO = 'U':
  *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
  *>
  *>             If IPIV(k) = IPIV(k-1) < 0, then rows and columns
  *>             k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
  *>             is a 2-by-2 diagonal block.
  *>
  *>          If UPLO = 'L':
  *>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  *>             interchanged and D(k,k) is a 1-by-1 diagonal block.
  *>
  *>             If IPIV(k) = IPIV(k+1) < 0, then rows and columns
  *>             k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
  *>             is a 2-by-2 diagonal block.
  *> \endverbatim
  *>
  *> \param[out] INFO
  *> \verbatim
  *>          INFO is INTEGER
  *>          = 0: successful exit
  *>          < 0: if INFO = -k, the k-th argument had an illegal value
  *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
  *>               has been completed, but the block diagonal matrix D is
  *>               exactly singular, and division by zero will occur if it
  *>               is used to solve a system of equations.
  *> \endverbatim
  *
  *  Authors:
  *  ========
  *
  *> \author Univ. of Tennessee
  *> \author Univ. of California Berkeley
  *> \author Univ. of Colorado Denver
  *> \author NAG Ltd.
  *
  *> \date December 2016
  *
  *> \ingroup complexSYcomputational
  *
  *> \par Further Details:
  *  =====================
  *>
  *> \verbatim
  *>
  *>  If UPLO = 'U', then A = U*D*U**T, where
  *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
  *>
  *>             (   I    v    0   )   k-s
  *>     U(k) =  (   0    I    0   )   s
  *>             (   0    0    I   )   n-k
  *>                k-s   s   n-k
  *>
  *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
  *>
  *>  If UPLO = 'L', then A = L*D*L**T, where
  *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
  *>
  *>             (   I    0     0   )  k-1
  *>     L(k) =  (   0    I     0   )  s
  *>             (   0    v     I   )  n-k-s+1
  *>                k-1   s  n-k-s+1
  *>
  *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
  *> \endverbatim
  *
  *> \par Contributors:
  *  ==================
  *>
  *> \verbatim
  *>
  *>  09-29-06 - patch from
  *>    Bobby Cheng, MathWorks
  *>
  *>    Replace l.209 and l.377
  *>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  *>    by
  *>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
  *>
  *>  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
  *>         Company
  *> \endverbatim
  *
  *  =====================================================================
        SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO )
  *
  *  -- LAPACK computational routine (version 3.7.0) --
  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  *     December 2016
  *
  *     .. Scalar Arguments ..
        CHARACTER          UPLO
        INTEGER            INFO, LDA, N
  *     ..
  *     .. Array Arguments ..
        INTEGER            IPIV( * )
        COMPLEX            A( LDA, * )
  *     ..
  *
  *  =====================================================================
  *
  *     .. Parameters ..
        REAL               ZERO, ONE
        PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
        REAL               EIGHT, SEVTEN
        PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
        COMPLEX            CONE
        PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
  *     ..
  *     .. Local Scalars ..
        LOGICAL            UPPER
        INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
        REAL               ABSAKK, ALPHA, COLMAX, ROWMAX
        COMPLEX            D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z
  *     ..
  *     .. External Functions ..
        LOGICAL            LSAME, SISNAN
        INTEGER            ICAMAX
        EXTERNAL           LSAME, ICAMAX, SISNAN
  *     ..
  *     .. External Subroutines ..
        EXTERNAL           CSCAL, CSWAP, CSYR, XERBLA
  *     ..
  *     .. Intrinsic Functions ..
        INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
  *     ..
  *     .. Statement Functions ..
        REAL               CABS1
  *     ..
  *     .. Statement Function definitions ..
        CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
  *     ..
  *     .. Executable Statements ..
  *
  *     Test the input parameters.
  *
        INFO = 0
        UPPER = LSAME( UPLO, 'U' )
        IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
           INFO = -1
        ELSE IF( N.LT.0 ) THEN
           INFO = -2
        ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
           INFO = -4
        END IF
        IF( INFO.NE.0 ) THEN
           CALL XERBLA( 'CSYTF2', -INFO )
           RETURN
        END IF
  *
  *     Initialize ALPHA for use in choosing pivot block size.
  *
        ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
  *
        IF( UPPER ) THEN
  *
  *        Factorize A as U*D*U**T using the upper triangle of A
  *
  *        K is the main loop index, decreasing from N to 1 in steps of
  *        1 or 2
  *
           K = N
     10    CONTINUE
  *
  *        If K < 1, exit from loop
  *
           IF( K.LT.1 )
       $      GO TO 70
           KSTEP = 1
  *
  *        Determine rows and columns to be interchanged and whether
  *        a 1-by-1 or 2-by-2 pivot block will be used
  *
           ABSAKK = CABS1( A( K, K ) )
  *
  *        IMAX is the row-index of the largest off-diagonal element in
  *        column K, and COLMAX is its absolute value.
  *        Determine both COLMAX and IMAX.
  *
           IF( K.GT.1 ) THEN
              IMAX = ICAMAX( K-1, A( 1, K ), 1 )
              COLMAX = CABS1( A( IMAX, K ) )
           ELSE
              COLMAX = ZERO
           END IF
  *
           IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN
  *
  *           Column K is zero or underflow, or contains a NaN:
  *           set INFO and continue
  *
              IF( INFO.EQ.0 )
       $         INFO = K
              KP = K
           ELSE
              IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  *
  *              no interchange, use 1-by-1 pivot block
  *
                 KP = K
              ELSE
  *
  *              JMAX is the column-index of the largest off-diagonal
  *              element in row IMAX, and ROWMAX is its absolute value
  *
                 JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
                 ROWMAX = CABS1( A( IMAX, JMAX ) )
                 IF( IMAX.GT.1 ) THEN
                    JMAX = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
                    ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
                 END IF
  *
                 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  *
  *                 no interchange, use 1-by-1 pivot block
  *
                    KP = K
                 ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
  *
  *                 interchange rows and columns K and IMAX, use 1-by-1
  *                 pivot block
  *
                    KP = IMAX
                 ELSE
  *
  *                 interchange rows and columns K-1 and IMAX, use 2-by-2
  *                 pivot block
  *
                    KP = IMAX
                    KSTEP = 2
                 END IF
              END IF
  *
              KK = K - KSTEP + 1
              IF( KP.NE.KK ) THEN
  *
  *              Interchange rows and columns KK and KP in the leading
  *              submatrix A(1:k,1:k)
  *
                 CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
                 CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
       $                     LDA )
                 T = A( KK, KK )
                 A( KK, KK ) = A( KP, KP )
                 A( KP, KP ) = T
                 IF( KSTEP.EQ.2 ) THEN
                    T = A( K-1, K )
                    A( K-1, K ) = A( KP, K )
                    A( KP, K ) = T
                 END IF
              END IF
  *
  *           Update the leading submatrix
  *
              IF( KSTEP.EQ.1 ) THEN
  *
  *              1-by-1 pivot block D(k): column k now holds
  *
  *              W(k) = U(k)*D(k)
  *
  *              where U(k) is the k-th column of U
  *
  *              Perform a rank-1 update of A(1:k-1,1:k-1) as
  *
  *              A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
  *
                 R1 = CONE / A( K, K )
                 CALL CSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
  *
  *              Store U(k) in column k
  *
                 CALL CSCAL( K-1, R1, A( 1, K ), 1 )
              ELSE
  *
  *              2-by-2 pivot block D(k): columns k and k-1 now hold
  *
  *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
  *
  *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
  *              of U
  *
  *              Perform a rank-2 update of A(1:k-2,1:k-2) as
  *
  *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
  *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
  *
                 IF( K.GT.2 ) THEN
  *
                    D12 = A( K-1, K )
                    D22 = A( K-1, K-1 ) / D12
                    D11 = A( K, K ) / D12
                    T = CONE / ( D11*D22-CONE )
                    D12 = T / D12
  *
                    DO 30 J = K - 2, 1, -1
                       WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
                       WK = D12*( D22*A( J, K )-A( J, K-1 ) )
                       DO 20 I = J, 1, -1
                          A( I, J ) = A( I, J ) - A( I, K )*WK -
       $                              A( I, K-1 )*WKM1
     20                CONTINUE
                       A( J, K ) = WK
                       A( J, K-1 ) = WKM1
     30             CONTINUE
  *
                 END IF
  *
              END IF
           END IF
  *
  *        Store details of the interchanges in IPIV
  *
           IF( KSTEP.EQ.1 ) THEN
              IPIV( K ) = KP
           ELSE
              IPIV( K ) = -KP
              IPIV( K-1 ) = -KP
           END IF
  *
  *        Decrease K and return to the start of the main loop
  *
           K = K - KSTEP
           GO TO 10
  *
        ELSE
  *
  *        Factorize A as L*D*L**T using the lower triangle of A
  *
  *        K is the main loop index, increasing from 1 to N in steps of
  *        1 or 2
  *
           K = 1
     40    CONTINUE
  *
  *        If K > N, exit from loop
  *
           IF( K.GT.N )
       $      GO TO 70
           KSTEP = 1
  *
  *        Determine rows and columns to be interchanged and whether
  *        a 1-by-1 or 2-by-2 pivot block will be used
  *
           ABSAKK = CABS1( A( K, K ) )
  *
  *        IMAX is the row-index of the largest off-diagonal element in
  *        column K, and COLMAX is its absolute value.
  *        Determine both COLMAX and IMAX.
  *
           IF( K.LT.N ) THEN
              IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
              COLMAX = CABS1( A( IMAX, K ) )
           ELSE
              COLMAX = ZERO
           END IF
  *
           IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN
  *
  *           Column K is zero or underflow, or contains a NaN:
  *           set INFO and continue
  *
              IF( INFO.EQ.0 )
       $         INFO = K
              KP = K
           ELSE
              IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  *
  *              no interchange, use 1-by-1 pivot block
  *
                 KP = K
              ELSE
  *
  *              JMAX is the column-index of the largest off-diagonal
  *              element in row IMAX, and ROWMAX is its absolute value
  *
                 JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
                 ROWMAX = CABS1( A( IMAX, JMAX ) )
                 IF( IMAX.LT.N ) THEN
                    JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
                    ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
                 END IF
  *
                 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  *
  *                 no interchange, use 1-by-1 pivot block
  *
                    KP = K
                 ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
  *
  *                 interchange rows and columns K and IMAX, use 1-by-1
  *                 pivot block
  *
                    KP = IMAX
                 ELSE
  *
  *                 interchange rows and columns K+1 and IMAX, use 2-by-2
  *                 pivot block
  *
                    KP = IMAX
                    KSTEP = 2
                 END IF
              END IF
  *
              KK = K + KSTEP - 1
              IF( KP.NE.KK ) THEN
  *
  *              Interchange rows and columns KK and KP in the trailing
  *              submatrix A(k:n,k:n)
  *
                 IF( KP.LT.N )
       $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
                 CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
       $                     LDA )
                 T = A( KK, KK )
                 A( KK, KK ) = A( KP, KP )
                 A( KP, KP ) = T
                 IF( KSTEP.EQ.2 ) THEN
                    T = A( K+1, K )
                    A( K+1, K ) = A( KP, K )
                    A( KP, K ) = T
                 END IF
              END IF
  *
  *           Update the trailing submatrix
  *
              IF( KSTEP.EQ.1 ) THEN
  *
  *              1-by-1 pivot block D(k): column k now holds
  *
  *              W(k) = L(k)*D(k)
  *
  *              where L(k) is the k-th column of L
  *
                 IF( K.LT.N ) THEN
  *
  *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
  *
  *                 A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
  *
                    R1 = CONE / A( K, K )
                    CALL CSYR( UPLO, N-K, -R1, A( K+1, K ), 1,
       $                       A( K+1, K+1 ), LDA )
  *
  *                 Store L(k) in column K
  *
                    CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
                 END IF
              ELSE
  *
  *              2-by-2 pivot block D(k)
  *
                 IF( K.LT.N-1 ) THEN
  *
  *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
  *
  *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
  *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
  *
  *                 where L(k) and L(k+1) are the k-th and (k+1)-th
  *                 columns of L
  *
                    D21 = A( K+1, K )
                    D11 = A( K+1, K+1 ) / D21
                    D22 = A( K, K ) / D21
                    T = CONE / ( D11*D22-CONE )
                    D21 = T / D21
  *
                    DO 60 J = K + 2, N
                       WK = D21*( D11*A( J, K )-A( J, K+1 ) )
                       WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
                       DO 50 I = J, N
                          A( I, J ) = A( I, J ) - A( I, K )*WK -
       $                              A( I, K+1 )*WKP1
     50                CONTINUE
                       A( J, K ) = WK
                       A( J, K+1 ) = WKP1
     60             CONTINUE
                 END IF
              END IF
           END IF
  *
  *        Store details of the interchanges in IPIV
  *
           IF( KSTEP.EQ.1 ) THEN
              IPIV( K ) = KP
           ELSE
              IPIV( K ) = -KP
              IPIV( K+1 ) = -KP
           END IF
  *
  *        Increase K and return to the start of the main loop
  *
           K = K + KSTEP
           GO TO 40
  *
        END IF
  *
     70 CONTINUE
        RETURN
  *
  *     End of CSYTF2
  *
        END
  