module Math.LinearAlgebra.Sparse.Algorithms.SolveLinear
(
solveLinear, solveLinSystems
)
where
import Data.Maybe
import Data.Monoid
import Data.IntMap as M hiding ((!))
import Math.LinearAlgebra.Sparse.Matrix
import Math.LinearAlgebra.Sparse.Vector
import Math.LinearAlgebra.Sparse.Algorithms.Staircase
import Math.LinearAlgebra.Sparse.Algorithms.Diagonal
solveDiagonal :: Integral α
=> (SparseMatrix α, SparseMatrix α, SparseMatrix α)
-> SparseVector α
-> Either String (SparseVector α)
solveDiagonal (d,t,u) b
| height d /= dim b = Left "width m /= dim b"
| isZeroVec b = Right $ zeroVec (width d)
| isZeroMx d = Left "left side is zero-matrix, while right side is not zero-vector"
| otherwise =
let (dd,a) = (mainDiag d, t ×· b)
(bad,sol) = solveDiag dd a
in if size (vec dd) < size (vec a) then Left "zeroes in lhs, while non-zeroes in rhs"
else if not (M.null bad) then Left $ "integral fraction error at lines "++(show (elems bad))
else Right $ (SV (width d) sol) ·× u
solveDiag :: Integral α => SparseVector α -> SparseVector α -> (IntMap Index, IntMap α)
solveDiag dd a = M.mapEitherWithKey solveOne (vec a)
where solveOne i r = let l = dd!i
(x,e) = r `divMod` l
in if (l == 0 && r /= 0) || e /= 0
then Left i else Right x
solveLinear :: Integral α =>SparseMatrix α -> SparseVector α -> Maybe (SparseVector α)
solveLinear m b = case solveDiagonal (toDiag m) b of
Left msg -> error msg
Right s -> Just s
solveLinSystems :: Integral α
=> SparseMatrix α
-> SparseVector (SparseVector α)
-> Maybe (SparseVector (SparseVector α))
solveLinSystems m bs = if ok then Just (SV (dim bs) sols) else Nothing
where (ok,sols) = M.mapAccum solve True (vec bs)
solve ok b = case solveLinear m b of
Nothing -> (False, emptyVec)
Just s -> (ok, s)