N=1000;
M=200;
K=256;
U=rand(N,M)+1i*rand(N,M);
V=rand(N,M)+1i*rand(N,M);
B=rand(N,K);
const_c=2;
const_d=0.5;
% no loops
tic
ux = U.'*B;
vx = V.'*B;
C_new = const_c*sum(abs(vx).^2,1);
D_new = const_d*sum(ux.*conj(vx),1);
toc
% with loops
tic
C=zeros(1,K);
D=zeros(1,K);
for m=1:M
ux=zeros(1,K);
vx=zeros(1,K);
for n=1:N
ux=ux+U(n,m)*B(n,:);
vx=vx+V(n,m)*B(n,:);
end
C=C+const_c*abs(vx).^2;
D=D+const_d*ux.*conj(vx);
end
toc
% results are not exactly the same, but the differences are ~1e-8 out of ~1e7:
[max(abs(C-C_new)) max(abs(D-D_new))]
% which is about as accurate as you can get for floating-point numbers around 1e7:
[max(eps(C)) max(eps(D))]
% plots to inspect the differences:
subplot(3,1,1)
plot(C,'--b','LineWidth',2)
hold on
plot(C_new,'-c')
title('C,C\_new')
subplot(3,1,2)
plot(real(D),'--r','LineWidth',2)
hold on
plot(real(D_new),'-m')
title('real(D,D\_new)')
subplot(3,1,3)
plot(imag(D),'--r','LineWidth',2)
hold on
plot(imag(D_new),'-m')
title('imag(D,D\_new)')
