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);