Optimize "xx = xi' * xi' and "S[k] += r[k,i] * xx" using BLAS
loope i and k reversed, some constants outside the loop
inner loop over i in logρ now a matrix multiplication
Big bad matmult for scatter matrix collection, use less memory in rnk()
Some cleaning up, scaling of the lower bound