def _gill_king(mat, beta, delta):
    """Backend function for the Gill-King algorithm."""
    size = mat.shape[0]

    # initialize d_vec and lowtri
    if scipy.sparse.issparse(mat):
        lowtri = scipy.sparse.eye(*mat.shape)
    else:
        lowtri = numpy.eye(size)

    d_vec = numpy.zeros(size, dtype=float)

    # there are no inner for loops, everything implemented with
    # vector operations for a reasonable level of efficiency
    for idx in range(size):
        if idx == 0:
            idz = []     # column index: all columns to left of diagonal
                        # d_vec(idz) doesn't work in case idz is empty
        else:
            idz = numpy.s_[:idx]

        djtemp = mat[idx, idx] - numpy.dot(
            lowtri[idx, idz], d_vec[idz]*lowtri[idx, idz].T)
        # C(idx, idx) in book

        if idx < size - 1:
            idy = numpy.s_[idx+1:size]
            # row index: all rows below diagonal
            ccol = mat[idy, idx] - numpy.dot(
                lowtri[idy, idz], d_vec[idz]*lowtri[idx, idz].T)
            # C(idy, idx) in book
            theta = abs(ccol).max()
            # guarantees d_vec(idx) not too small and lowtri(idy, idx) not too
            # big in sufficiently positive definite case, d_vec(idx) = djtemp
            d_vec[idx] = max(abs(djtemp), (theta/beta)**2, delta)
            lowtri[idy, idx] = ccol/d_vec[idx]

        else:
            d_vec[idx] = max(abs(djtemp), delta)

    # convert to usual output format: replace lowtri by lowtri*sqrt(D) and
    # transpose
    for idx in range(size):
        lowtri[:, idx] = lowtri[:, idx]*numpy.sqrt(d_vec[idx])
        # lowtri = lowtri*diag(sqrt(d_vec)) bad in sparse case

    return lowtri
Backend function for the Gill-King algorithm.