function Z = array3mult(X,A,B,C)
% Z = array3mult(X,A,B,C)
%
% Z is a 3-dimemsional array such that
%   Z(:) = kron(C,kron(B,A))*X(:)

% Per Christian Hansen, IMM, August 21, 2005.

[m,n,p] = size(X);

if size(A,2)~=m, error('Incorrect size of A'), end
if size(B,2)~=n, error('Incorrect size of B'), end
if size(C,2)~=p, error('Incorrect size of C'), end

% Multiplication along first and second dimensions.
Z = zeros(m,n,p);
for i=1:p, Z(:,:,i) = A*X(:,:,i)*B'; end

% Multiplication along third dimension.
Z = shiftdim(Z,2);
for i=1:n, Z(:,:,i) = C*Z(:,:,i); end
Z = shiftdim(Z,1);