Вы можете использовать outer с размерами матрицы.
outer
n <- 10 p <- 5 A <- matrix( rnorm(n*p), n, p ) f <- function(x,y) sqrt(sum((x-y)^2)) B <- outer( 1:n, 1:n, Vectorize( function(i,j) f(A[i,], A[j,]) ) )