You could take this approach:
A = your matrix
AA = A*A'; % all of the dot products between rows (more than you need, actually)
d = sqrt(diag(AA)); % the norms of the rows
AA = AA ./ d; AA = AA ./ d'; % normalize everything, or AA = AA ./ (d*d');
The result has the cosine similarities between all of the rows. Pick off what you need. E.g., AA(1,2) has the cosine similarity between rows 1 and 2. AA(5,8) has the cosine similarity between rows 5 and 8. Etc.
This does extra computations that you don't need, of course. But by using functions like matrix multiply that are multi-threaded I am guessing it may run faster than custom code that tries to avoid the unnecessary computations. The caveat is if the row size is too large then the memory usage may bog you down.
