'A function used by CVXOPT\'s nonlinear solver.
Based on f(x,y,z) = y * log(y) + x - y * log(z).
Parameters
vars_: A cvxopt dense matrix with values for (x,y,z).
scaling: A scaling for the Hessian.
Returns
_solver_hook() returns the constraint size and a feasible point.
_solver_hook(x) returns the function value and gradient at x.
_solver_hook(x, z) returns the function value, gradient,
and (z scaled) Hessian at x.'

def _solver_hook(self, vars_=None, scaling=None):

    import cvxopt
    entries = (self.size[0] * self.size[1])
    if (vars_ is None):
        x_init = (entries * [0.0])
        y_init = (entries * [0.5])
        z_init = (entries * [1.0])
        return (self.size[0], cvxopt.matrix(((x_init + y_init) + z_init)))
    x = vars_[0:entries]
    y = vars_[entries:(2 * entries)]
    z = vars_[(2 * entries):]
    if ((min(y) <= 0.0) or (min(z) <= 0.0)):
        return None
    f = cvxopt.matrix(0.0, (entries, 1))
    for i in range(entries):
        f[i] = ((x[i] - (y[i] * math.log(z[i]))) + (y[i] * math.log(y[i])))
    Df = cvxopt.matrix(0.0, (entries, (3 * entries)))
    for i in range(entries):
        Df[(i, i)] = 1.0
        Df[(i, (entries + i))] = ((math.log(y[i]) - math.log(z[i])) + 1.0)
        Df[(i, ((2 * entries) + i))] = ((- y[i]) / z[i])
    if (scaling is None):
        return (f, Df)
    big_H = cvxopt.spmatrix(0, [], [], size=((3 * entries), (3 * entries)))
    for i in range(entries):
        H = cvxopt.matrix([[0.0, 0.0, 0.0], [0.0, (1.0 / y[i]), ((-1.0) / z[i])], [0.0, ((-1.0) / z[i]), (y[i] / (z[i] ** 2))]])
        big_H[i:(3 * entries):entries, i:(3 * entries):entries] = (scaling[i] * H)
    return (f, Df, big_H)
