function D = dmatrix(P)

N = size(P,1);

d = zeros(N);
for i = 1:(N-1)
    X = ones((N-i),1)*P(i,:);
    d((i+1):N,i) = sum((X - P((i+1):N,:)).^2,2);
end

D = sqrt(d + d');