 AA=[
   0.3210             0.1743 - 0.0657i  -0.3102 + 0.0027i
   0.1743 + 0.0657i   0.2081            -0.0934 - 0.1455i
  -0.3102 - 0.0027i  -0.0934 + 0.1455i   0.4709          ];
  AA=AA/trace(AA);
  [ua,da]=eig(AA);
  da=real(da);
  eiga=sort(diag(da),'descend');
 BB=[
   0.1304            -0.0548 - 0.0460i  -0.0564 - 0.1229i
  -0.0548 + 0.0460i   0.5758             0.0976 - 0.1704i
  -0.0564 + 0.1229i   0.0976 + 0.1704i   0.2937          ];
  BB=BB/trace(BB);
  [ub,db]=eig(BB);
  db=real(db);
  eigb=sort(diag(db),'descend');
  [eiga eigb]
  BB=ua*db*ua';   % now commuting
  B{1}=BB;
  A{1}=AA;
  %norm(AA*BB-BB*AA)
%%%%take convex comb. of n random permutation matrices
%%%In=eye(n);
%%%xxx=rand(n,1)+1;
%%%xxx=xxx/sum(xxx);
%%%DS=zeros(n);
%%%for iii=1:n,
%%%   DS=DS+xxx(iii)*In(randperm(n),:);
%%%end
%%%eiga=rand(n,1)+1;
%%%eiga=eiga/sum(eiga);
%%%eigb=DS*eiga;    % eiga larger majorization, eigb < eiga
%%%[q,~]=qr(randn(n));
%%%AA=q*diag(eiga)*q';
%%%BB=q*diag(eigb)*q';
%%%B{1}=BB;
%%%A{1}=AA;
