def constraint(u, H, gamma, D, C, tol, verbose = 0):
    """Check the constraint violation for given u, fixed N

    Parameters
    ----------
    u : np.array of shape(None,)
        Beamlet radiation distribution 
    H : Sparse array of shape(None, None)
        Block-diagonal OAR dose-deposition matrix block_diag(H_1, ... , H_M) 
        where H_i is the sparse dose deposition matrix for i-th modality (for a given OAR)
    gamma : np.array of shape (None,)
        Linear block coefficient array of the OAR BE 
        (with N_m included as coefficients, vector gamma_tilde in the paper)
    D : np.array of shape (None,)
        Quadratic block-diagonal coefficient array D of the OAR BE
        (with N_m included as coefficients, D as it is in the paper, currently the notation is not ideal)
    C : float
        Constraint to satisfy
    tol : float
        Tolerance to check the soft constraint violation
    verbose : 0 or 1
        0 for no printing
        1 for printing
   
    Returns
    -------
    tuple (bool, float, float, float)
        (soft_constr_satisfied, distance, constr_val at u, C)
        
    Notes
    -----
    Changes: np.diag(some_array)@another_array changed to some_array*another_array
             changed G notation to D notation as per the paper
             changed B notation to H notation as per the paper
             changed the order of the arguments to more logical
             changed to the possibly sparse H interface, @ -> H.dot
        
    """
    x = H.dot(u)
    constr_val = gamma.T.dot(x) + x.T.dot(D*x)
    hard = constr_val <= C
    soft = constr_val <= (1+tol)*C
    if bool(verbose):
        print('   Hard constraint satisfied: {}'.format(hard))
        print('   Relaxed constraint within {}% tol satisfied: {}'.format(tol*100, soft))
    return soft, 100*(constr_val - C)/C, constr_val, C