How to implement a constraint in a p-hub location problem?

I try to implement two constraints from a p-hub median location problem.The constraints are :

Z_ij^km ≤ x_m for all i,j,k,m

Z_ij^km ≤ x_k for all i,j,k,m

Z is a 4-D (i*j*k*m) array and x a vector. Z and x are my unknown. In my case, i=j=k=m=s=20 (but, if it's possible, it should go up to 100)

I tried this code but, first, I'm not sure if it's exact and it is incredibly slow.

clearer1 = zeros(size(objX));
clearer12 = clearer1(:);
clearer2 = zeros(size(objZ));
clearer22 = clearer2(:);

% ----Inequality constraints----

Aineq = spalloc(2*s*s*s*s,s2,4*s*s*s*s);
bineq = spalloc(2*s*s*s*s,1,0);

counter = 1;

% On x_m
for i=1:s
    for j=1:s
        for k=1:s
            for m=1:s
                Ztemp = clearer2;
                Ztemp(i,j,k,m) = 1;
                Xtemp = clearer1;
                Xtemp(m) = -1;
                Temp = sparse([Xtemp(:);Ztemp(:)]);
                Aineq(counter,:) = Temp';
                counter = counter+1;
            end
        end
    end
end

% On x_k
for i=1:s
    for j=1:s
        for k=1:s
            for m=1:s
                Ztemp = clearer2;
                Ztemp(i,j,k,m) = 1;
                Xtemp = clearer1;
                Xtemp(k) = -1;
                Temp = sparse([Xtemp(:);Ztemp(:)]);
                Aineq(counter,:) = Temp';
                counter = counter+1;
            end
        end
    end
end

How can I write correctly my constraints? How can I improve the speed of my function?