function S = skron(A,B);
% symmetrized Kronecker product
%
% skron is a matrix representation for the linear operator from S_n to S_n:
%        K    ->    1/2 (BKA' + AKB')
% so that
%      skron(A,B) svec(K) = 1/2 svec(BKA' + AKB').
% where svec maps a symmetric matrix to a vector of length n(n+1)/2,
% multiplying off diagonal components by sqrt(2), so svec(K)'svec(L) = tr KL.
% This last property implies that svec is a unitary operator.

% kron is a matrix representation for the linear operator from M_n to M_n:
%        K    ->    BKA'
% so that 
%      kron(A,B) vec(K) = vec(BKA')
% and therefore
%      (kron(A,B) + kron(B,A)) vec(K) = vec(BKA' + AKB').
% where vec maps a square matrix to a vector of length n^2.
% We make use of the identity
%
%      skron(A,B) = svec mat 1/2 (kron(A,B) + kron(B,A)) vec smat
%                 = svec mat 1/2 (kron(A,B) + kron(B,A)) mat^* svec^*
%                 = svec sym mat 1/2 (kron(A,B) + kron(B,A)) mat^* sym^* svec^*
% where mat is the inverse of vec (also the adjoint of vec) and
%      smat is the inverse of svec (also the adjoint of svec).
%       sym symmetrizes a matrix (sym^* embeds a symmetric matrix in M_n)
%
% Thus to compute skron(A,B), we start with 1/2 (kron(A,B) + kron(B,A))
% (which is provided by the Matlab built-in function kron), then apply
% the operator "svec sym mat" to each of its columns, and then apply
% the same operator to the transposes of each of the resulting rows, and 
% transpose the result. 
% Copyright J.V. Burke, A.S. Lewis, M.L. Overton, 2001.

n = length(A);
K = 0.5*(kron(A,B) + kron(B,A));
for k = 1: n^2
   col = K(:,k);
   M = reshape(col,n,n);
   M = 0.5*(M + M');
   L(:,k) = svec(M,n);   
end
for k = 1: n*(n+1)/2
   row = L(k,:);
   M = reshape(row',n,n);
   M = 0.5*(M + M');
   S(k,:) = (svec(M,n))';   
end
