'Solve quadratic subproblem'

def solve(self, tr_radius):

    (lambda_current, lambda_lb, lambda_ub) = self._initial_values(tr_radius)
    n = self.dimension
    hits_boundary = True
    already_factorized = False
    self.niter = 0
    while True:
        if already_factorized:
            already_factorized = False
        else:
            H = (self.hess + (lambda_current * np.eye(n)))
            (U, info) = self.cholesky(H, lower=False, overwrite_a=False, clean=True)
        self.niter += 1
        if ((info == 0) and (self.jac_mag > self.CLOSE_TO_ZERO)):
            p = cho_solve((U, False), (- self.jac))
            p_norm = norm(p)
            if ((p_norm <= tr_radius) and (lambda_current == 0)):
                hits_boundary = False
                break
            w = solve_triangular(U, p, trans='T')
            w_norm = norm(w)
            delta_lambda = ((((p_norm / w_norm) ** 2) * (p_norm - tr_radius)) / tr_radius)
            lambda_new = (lambda_current + delta_lambda)
            if (p_norm < tr_radius):
                (s_min, z_min) = estimate_smallest_singular_value(U)
                (ta, tb) = self.get_boundaries_intersections(p, z_min, tr_radius)
                step_len = min([ta, tb], key=abs)
                quadratic_term = np.dot(p, np.dot(H, p))
                relative_error = (((step_len ** 2) * (s_min ** 2)) / (quadratic_term + (lambda_current * (tr_radius ** 2))))
                if (relative_error <= self.k_hard):
                    p += (step_len * z_min)
                    break
                lambda_ub = lambda_current
                lambda_lb = max(lambda_lb, (lambda_current - (s_min ** 2)))
                H = (self.hess + (lambda_new * np.eye(n)))
                (c, info) = self.cholesky(H, lower=False, overwrite_a=False, clean=True)
                if (info == 0):
                    lambda_current = lambda_new
                    already_factorized = True
                else:
                    lambda_lb = max(lambda_lb, lambda_new)
                    lambda_current = max(np.sqrt((lambda_lb * lambda_ub)), (lambda_lb + (self.UPDATE_COEFF * (lambda_ub - lambda_lb))))
            else:
                relative_error = (abs((p_norm - tr_radius)) / tr_radius)
                if (relative_error <= self.k_easy):
                    break
                lambda_lb = lambda_current
                lambda_current = lambda_new
        elif ((info == 0) and (self.jac_mag <= self.CLOSE_TO_ZERO)):
            if (lambda_current == 0):
                p = np.zeros(n)
                hits_boundary = False
                break
            (s_min, z_min) = estimate_smallest_singular_value(U)
            step_len = tr_radius
            if (((step_len ** 2) * (s_min ** 2)) <= ((self.k_hard * lambda_current) * (tr_radius ** 2))):
                p = (step_len * z_min)
                break
            lambda_ub = lambda_current
            lambda_lb = max(lambda_lb, (lambda_current - (s_min ** 2)))
            lambda_current = max(np.sqrt((lambda_lb * lambda_ub)), (lambda_lb + (self.UPDATE_COEFF * (lambda_ub - lambda_lb))))
        else:
            (delta, v) = singular_leading_submatrix(H, U, info)
            v_norm = norm(v)
            lambda_lb = max(lambda_lb, (lambda_current + (delta / (v_norm ** 2))))
            lambda_current = max(np.sqrt((lambda_lb * lambda_ub)), (lambda_lb + (self.UPDATE_COEFF * (lambda_ub - lambda_lb))))
    self.lambda_lb = lambda_lb
    self.lambda_current = lambda_current
    self.previous_tr_radius = tr_radius
    return (p, hits_boundary)
