Due to Fortran's memory model, its performance could potentially be drastically improved by altering the arrays such that the program is accessing consecutive columns rather than rows.

from: r = r + (X(i,k) - X(j,k)) * (X(i,k) - X(j,k))

to: r = r + (X(k,i) - X(k,j)) * (X(k,i) - X(k,j))