"""
Higams' ldlt_symm translated from matlab. Performs  a so called
 BKK  Bounded Bunch Kaufman factorization of A0, that means 
 ||L|| is bounded and bounded away from zero.
"""
function  ldlt_symm(A0 :: Array{Float64,2}, piv :: Char='r')

    #LDLT_SYMM  Block LDL^T factorization for a symmetric indefinite matrix.
    #     Given a Hermitian matrix A,
    #     [L, D, P, RHO, NCOMP] = LDLT_SYMM(A, PIV) computes a permutation P,
    #     a unit lower triangular L, and a real block diagonal D
    #     with 1x1 and 2x2 diagonal blocks, such that  P*A*P' = L*D*L'.
    #     PIV controls the pivoting strategy:
    #       PIV = 'p': partial pivoting (Bunch and Kaufman),
    #       PIV = 'r': rook pivoting (Ashcraft, Grimes and Lewis).
    #     The default is partial pivoting.
    #     RHO is the growth factor.
    #     For PIV = 'r' only, NCOMP is the total number of comparisons.
    
    #     References:
    #     J. R. Bunch and L. Kaufman, Some stable methods for calculating
    #        inertia and solving symmetric linear systems, Math. Comp.,
    #        31(137):163-179, 1977.
    #     C. Ashcraft, R. G. Grimes and J. G. Lewis, Accurate symmetric
    #        indefinite linear equation solvers. SIAM J. Matrix Anal. Appl.,
    #        20(2):513-561, 1998.
    #     N. J. Higham, Accuracy and Stability of Numerical Algorithms,
    #        Second edition, Society for Industrial and Applied Mathematics,
    #        Philadelphia, PA, 2002; chap. 11.
    
    #    This routine does not exploit symmetry and is not designed to be
    #     efficient.
    #
    #    Adapted from N. Higham Matlab code.
    #    JPD april 23 2015


    # Since array contents are mutable and modified,
    # copy the input matrix to keep it unchanged outside the
    # scope of the function

    A = copy(A0)

    # minimal checks for conforming inputs
    isequal(triu(A)',tril(A)) || error("Must supply Hermitian matrix.")
    piv in ['p','r'] || error("Pivoting must be 'p' or 'r'.")

    n, = size(A)

    k = 1
    D = eye(n,n); 
    L = eye(n,n);  
    if n == 1   D = A; end
    
    pp = collect(1:n)
    
    maxA = norm(A, Inf)
    ρ = maxA;
    
    ncomp = 0;
    s=1
    
    α = (1 + sqrt(17))/8
    while k < n
        (λ, vr) = findmax( abs.(A[k+1:n,k]) );
        r = vr[1] + k;
        if λ > 0
            swap = false;
            if abs(A[k,k]) >= α*λ
                s = 1;
            else 
                if piv == 'p'
                    temp = A[k:n,r]; temp[r-k+1] = 0
                    σ = norm(temp, Inf)
                    if α*λ^2 <= abs(A[k,k])*σ
                        s = 1;
                    elseif abs(A[r,r]) >= α*σ
                        swap = true;
                        m1 = k; m2 = r;
                        s = 1;
                    else
                        swap = true;
                        m1 = k+1; m2 = r;
                        s = 2;
                    end
                    if swap
                        A[[m1,m2],:] = A[[m2,m1],:]
                        L[[m1,m2],:] = L[[m2,m1],:]
                        A[:,[m1,m2]] = A[:,[m2,m1]]
                        L[:,[m1,m2]] = L[:,[m2,m1]]
                        
                        pp[[m1,m2]] = pp[[m2,m1]]
                    end
                elseif piv == 'r'
                    j = k;
                    pivot = false;
                    λ_j = λ;
                    while ~pivot
                        (temp1,vr) = findmax( abs.(A[k:n,j]) );
                        ncomp = ncomp + n-k;
                        r = vr[1] + k - 1;
                        temp = A[k:n,r]; temp[r-k+1] = 0.0;
                        λᵣ = max(maximum(temp), -minimum(temp))
                        ncomp = ncomp + n-k;
                        if α*λᵣ <= abs(A[r,r])
                            pivot = true;
                            s = 1;
                            A[k,:], A[r,:] = A[r,:], A[k,:]
                            L[k,:], L[r,:] = L[r,:], L[k,:]
                            A[:,k], A[:,r] = A[:,r], A[:,k]
                            L[:,k], L[:,r] = L[:,r], L[:,k]
                            pp[k], pp[r]  = pp[r], pp[k]
                        elseif λ_j == λᵣ
                            pivot = true;
                            s = 2;
                            A[k,:], A[j,:] = A[j,:], A[k,:]
                            L[k,:], L[j,:] = L[j,:], L[k,:]
                            A[:,k], A[:,j] = A[:,j], A[:,k]
                            L[:,k], L[:,j] = L[:,j], L[:,k]
                            pp[k], pp[j]  = pp[j], pp[k]
                            k1 = k+1;
                            A[k1,:], A[r,:] = A[r,:], A[k1,:]
                            L[k1,:], L[r,:] = L[r,:], L[k1,:]
                            A[:,k1], A[:,r] = A[:,r], A[:,k1]
                            L[:,k1], L[:,r] = L[:,r], L[:,k1]
                            pp[k1], pp[r]  = pp[r], pp[k1]
                        else
                            j = r;
                            s = 1
                            λ_j = λᵣ;
                        end
                    end
                end
            end
            
            if s == 1
                
                D[k,k] = A[k,k];
                A[k+1:n,k] = A[k+1:n,k]/A[k,k];
                L[k+1:n,k] = A[k+1:n,k];
                i = k+1:n;
                A[i,i] = A[i,i] - A[i,k:k] * A[k:k,i];
                A[i,i] = 0.5 * (A[i,i] + A[i,i]');
                
            elseif s == 2
                
                E = A[k:k+1,k:k+1];
                D[k:k+1,k:k+1] = E;
                i= k+2:n;
                C = A[i,k:k+1];
                temp = C/E;
                L[i,k:k+1] = temp;
                A[i,k+2:n] = A[i,k+2:n] - temp*C';
                A[i,i] = 0.5 * (A[i,i] + A[i,i]');                         
            end
            
            # Ensure diagonal real (see LINPACK User's Guide, p. 5.17).
            for i=k+s:n
                A[i,i] = real(A[i,i]);
            end
            
            if  k+s <= n
                val, = findmax(abs.(A[k+s:n,k+s:n]))
                ρ = max(ρ, val );
            end
            
        else  # Nothing to do.           
            s = 1;
            D[k,k] = A[k,k];            
        end
        
        k = k + s;
        
        if k == n
            D[n,n] = A[n,n];
            break
        end
        
    end
    #P=eye(n)
    #P = P[pp,:]
    return L, D, pp, ρ, ncomp
end
