function y = afun(x, B, C, G, H, m, n)
% Compute the product
%   y = A x
% where
%   B = (B)_{ij}, i,j = 1 : m
%   C = (C)_{ij}
%   A = (A)_{ij}
%   (A)_{ij} = (B)_{ij} \kron (C)_{ij}
% Input:
%   x : vector, mn^2
%   B : matrix, mn-by-mn
%   C : matrix, mn-by-mn
%   G, H : ignored
%   m : scalar, >= 1
%   n : scalar, >= 1
% Output:
%   y : vector, mn^2

if nargin < 7
   error(sprintf('afun error: only %d arguments, need 7', nargin))
end

n2 = n*n;
mn2 = m * n2;
mn = m * n;

% Check sizes.

[mB, nB] = size(B);
[mC, nC] = size(C);
[mx, nx] = size(x);

if (mB ~= mn) | (nB ~= mn)
   error('afun error: B must be mn-by-mn')
end
if (mC ~= mn) | (nC ~= mn)
   error('afun error: C must be mn-by-mn')
end
if (mx ~= mn2) | (nx ~= 1)
   error('afun error: x must be mn^2-by-1')
end

% Do multiplication.

y = zeros(mn2,1);

for i = 1 : m
   % Compute (y)_i = y( (i-1)*n^2 + (1:n^2) ).
   ii1 = (i-1)*n + (1:n);
   ii2 = (i-1)*n2 + (1:n2);
   yy = zeros(n2,1);
   for j = 1 : m
      jj1 = (j-1)*n + (1:n);
      jj2 = (j-1)*n2 + (1:n2);
      yy = yy + KPmultV(B(ii1,jj1), C(ii1,jj1), x(jj2));
   end
   y(ii2) = yy;
end
