module Math.LinearAlgebra.Sparse.Algorithms.Diagonal
( 
isDiag, toDiag
)
where

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 

-- | Checks if matrix has diagonal form
isDiag :: SparseMatrix α -> Bool
isDiag m = M.foldlWithKey' (\true i row -> true && (keys row == [i])) True (mx m)

-- | Transforms matrix to diagonal form and returns also two protocol matrices:
--
-- >>> let (d,t,u) = toDiag m  in t × m × (trans u) == d
-- True
--
--   @t@ stores rows transformations and @u@ — columns transformations
toDiag :: Integral α => SparseMatrix α -> (SparseMatrix α, SparseMatrix α, SparseMatrix α)
toDiag m = toDiag' m (idMx (height m))
                     (idMx (width  m))

toDiag' m t u | isZeroMx m =  (m,t,u)
              | height m /= height t = error "height(mM) /= height(mT)"
              | width m /= width u   = error "width(mM) /= width(mU)"
              | otherwise = 
              let (s,t') = staircase' m t
              in  dm True s t' u

-- Here  s  is a staircase matrix.
-- If it is not diagonal, then transp(s) it brought to staircase
-- form - this corresponds to the column elementary 
-- transformations of s -  and so on,  until the diagonal matrix 
-- is obtained (even number of `transp' to be applied).
dm even s t u = case (even, isDiag s, trans s) of
    (True , True , _ ) -> (s , t, u)
    (False, True , s') -> (s', t, u)
    (True , False, s') -> dm False s'' t  u' where (s'',u') = staircase' s' u
    (False, False, s') -> dm True  s'' t' u  where (s'',t') = staircase' s' t