% bprod.m for computing the biproduct  A \odot B
% Usage:     C = bprod(A,B)
%
% e.g. A=diag(1:4); I=eye(4); C=bprod(A,I), eig(A), eig(2*C)
%  or  A=round(rand(4)*9+1); A=tril(A,1); A=tril(A',1);A=sparse(A);
%      I=speye(4); C=bprod(A,I), eigs(A),eigs(2*C)
%      hold off,spy(A), hold on, spy(C,'dr',13); hold off
%% K Chen (2004)

function C = bprod(A,B);
if nargin<2, help bprod, return, end
if size(A)~=size(B), help bprod, disp('CHECK sizes!'); return, end

n=size(A,1); C=[];

if ~issparse(A) & ~issparse(B)  % ---------------- Non-sparse case 

for p=2:n
for r=2:n
  for q=1:p-1
  for s=1:r-1
   u = (p-2)*(p-1)/2 + q; v = (r-2)*(r-1)/2 + s; 
   C(u,v) = (A(p,r)*B(q,s)-A(p,s)*B(q,r) + ...
             B(p,r)*A(q,s)-B(p,s)*A(q,r) ) / 2;
  end
  end
end
end

else %----------------------------------------------- Sparse below

   % [ia,ja,sa] = find(A); [ib,jb,sb] = find(B); Not optimal below!!!
   C=sparse(n,n);
for p=2:n
for r=2:n
  for q=1:p-1
  for s=1:r-1
   u = (p-2)*(p-1)/2 + q; v = (r-2)*(r-1)/2 + s; 
   C(u,v) = (A(p,r)*B(q,s)-A(p,s)*B(q,r) + ...
             B(p,r)*A(q,s)-B(p,s)*A(q,r) ) / 2;
   if abs(C(u,v))<eps,C(u,v)=sparse(1,1);end
  end
  end
end
end

end  %----------------------------------------------- |=|
