haskell-numerics / hmatrix

Linear algebra and numerical computation
381 stars 104 forks source link

`triDiagSolve` mutates its arguments #282

Closed menelaos closed 5 years ago

menelaos commented 6 years ago

triDiagSolve mutates its first and third argument (corresponding to the lower and upper diagonal). Here is a minimal(-ish) working example:

import Numeric.LinearAlgebra

-- Solve kx = b via `linearSolve`
k = matrix 3 
  [  2, -1,  0
  , -1,  2, -1
  ,  0, -1,  2
  ]

b = col [10, 10, 10]

correctSolution = linearSolve k b

-- Solve the same equation with `triDiagSolve`
minusOnes = vector [-1, -1]
twos      = vector [2, 2, 2]

wrongSolution = triDiagSolve minusOnes twos minusOnes b

main = do
  putStrLn $ "This is the (correct) solution found by `linearSolve`:"
  print correctSolution
  putStrLn $ "This is the (incorrect) solution found by `triDiagSolve`:"
  print wrongSolution
  putStrLn $ "This should be a vector of -1s:"
  print minusOnes
  -- [-0.5,-0.5714285714285714]

In the example above, triDiagSolve produces an incorrect result because the same vector (minusOnes) is used for both the upper and lower diagonal and gets mutated along the way.

With k and b given as above, only the first argument is mutated, but for larger matrices, the third argument is affected as well.

idontgetoutmuch commented 6 years ago

Thanks for the report. It won't get fixed in the short term by me but PRs are always welcome :-)