dense-linear-algebra : Weird memory and runtime behavior from `generateSym`

[Magalame]

The generateSym function is defined as:

generateSym :: Int -> (Int -> Int -> Double) -> Matrix
generateSym n f = runST $ do
  m <- unsafeNew n n
  for 0 n $ \r -> do
    unsafeWrite m r r (f r r)
    for (r+1) n $ \c -> do
      let x = f r c
      unsafeWrite m r c x
      unsafeWrite m c r x
  unsafeFreeze m

Running it with n=100, I noted we can note that the function allocates ~ 160 000 bytes of memory, which is around twice what we would expect when allocating one Matrix. This allocation seems to be related to the dependence on c of x. If we change f r c to f r r, the allocation drops 80 000 bytes, and the runtime is divided by two.

[ocramz]

If we change f r c to f r r

Doesn't this change the behaviour of the function as well? i.e. it would only read indices from the diagonal. I think it would be better to introduce sparse formats, such as "lower triangular" matrices etc. and have specialized implementations for those. But this would also significantly increase the scope of this library.

[Magalame]

Doesn't this change the behaviour of the function as well?

It does, it was just to point out that there is something weird happening with the dependence on c, because there seems to be no obvious reason why f r c would lead to twice the allocation induced by f r r.

But this would also significantly increase the scope of this library.

At some point we'll probably need these features in the Haskell ecosystem anyway :)