20120702 Libing wang recommend:
[nSmp nFea]=size(fea);
Dis = pdist2(fea,fea);
------------------------------------------【以下可不看】 ------------------------------------------- 
%能不用循环尽量不用循环
[nSmp nFea]=size(fea);
normfea = sum(fea .^2,2);
Dis = -2*fea*fea'+repmat(normfea,1,nSmp)+repmat(normfea',nSmp,1);
Dis(Dis<=0)  = 0;%有次程序运行，发现对角线上有很小的非零的值，应该是计算精度的问题
Dis = sqrt(Dis);