function C = mat_mul_mat_3(A, B)
% C = mat_mul_mat(A, B):           C = A*B
% multiplication with N x 3 matrix B
C = zeros(size(B));
if ndims(A) == 3
for i=1:3
    for k=1:3
        C(:,i) = C(:,i) + A(:,i,k).*B(:,k);
    end
end
else
for i=1:3
    for k=1:3
        C(:,i) = C(:,i) + A(i,k).*B(:,k);
    end
end
end