This function computes pairwise distance between two sample sets and produce a matrix of square of Euclidean or Mahalanobis distances. The code is fully optimized by vectorization. Therefore it is much faster than the built-in function pdist.
When two matrices A and B are provided as input, this function computes the square Euclidean distances between them. If an extra positive definite matrix M is provided, it computes Mahalanobis distances.

If only one matrix A is provided, the function computes pairwise square Euclidean distances between vectors in A. In this case, it is equivalent to the square of pdist function in matlab statistics toolbox but much faster.

Sample code:
d=1000;n1=5000;n2=6000;
A=rand(d,n1);B=rand(d,n2);
M=rand(d,d);M=M*M'+eye(d);
D1=sqdist(A,B);
D2=sqdist(A);
D3=sqdist(A,B,M);

Detail explanation can be found in following blog post:
http://statinfer.wordpress.com/2011/11/14/efficient-matlab-i-pairwise-distances/

This function is now a part of the PRML toolbox (http://www.mathworks.com/matlabcentral/fileexchange/55826-pattern-recognition-and-machine-learning-toolbox).