The problem is a standard nonlinear nonconvex problem, so any solver for this problem class is suitable to solve the problem.
As an example, the following code implements the problem in the MATLAB Toolbox YALMIP (disclaimer, developed by me) and solves the problem using the local nonlinear solver ipopt.
N = 3;
K = 3;
h = rand(N,K);
R = rand(K,1);
sigma = 1;
beta = 1;
Pmax = 10;
Constraints = ;
P = sdpvar(N,K,'full');
C = sdpvar(N,K,'full');
Constraints = [P>=0, C>=0, sum(sum(P)) <= Pmax,sum(C,1)<=1]
for k = 1:K
Constraints = [Constraints, t<=(beta/(N*R(k)))*sum(C(:,k).*log2(1+P(:,k).*h(:,k)/sigma))];
ipopt is a local solver, so all you can hope for here is a locally optimal solution. However, it looks as if your problem actually is rather easy, as the solution often seems to be globally optimal. You can see this by solving the problem globally using the global solver available in YALMIP. It implements a b&b strategy based on linear relaxations for lower bounds and nonlinear local solvers for upper bounds. Of course, it is only applicable to small problems (unless you have a lot of time to spare...)
Finally, I would like to add that the structure of your problem actually allows you to eliminate the c-variables. It can be seen that the optimal choice of c is to have all elements in each column equal to zero, execpt the element which corresponds to the largest term in the vector of logarithmic terms. Hence, conceptually, you can solve
Constraints = [P>=0, sum(sum(P))<= Pmax]
for k = 1:K
Constraints = [Constraints, t<=(beta/(N*R(k)))*max(log2(1+P(:,k).*h(:,k)/sigma))];
This leads to a mixed-integer nonlinear nonconvex program when you model it in YALMIP (binary variables are introduced to model the nonconvex use of the max operator), and it seems to be a less efficient approach than the straightforward nonlinear approach.