module Matrix.Pivot (split1, permute, permuteLL, vswap, mswap) where 
--TODO Find better home for split1 and maybe others   


import Data.Array
import Data.List

--
-- Trivial Pivoting Routines
-- 
mswap :: [Int] -> Array (Int,Int) Double -- ^ A
   -> Array (Int,Int) Double -- ^ LU(A)

mswap n a = a'
    where a' = array bnds [ ((i,j), luij i j) | (i,j) <- range bnds ]
	  luij i j = a!(ie i,j)
	  bnds = bounds a
	  ie i = n!!(i - 1)

vswap :: [Int] -> Array Int Double -- ^ A
   -> Array Int Double -- ^ LU(A)

vswap n a = a'
    where a' = array bnds [ (j, luij j) | j <- range bnds ]
	  luij j = a!(ie j)
	  bnds = bounds a
	  ie j = n!!(j - 1)   
	  


rows1 :: Int -> [Double] -> [Int]
rows1 _ [] = []
rows1 n (x:xs) | abs (x) > 0 = (n + 1) : rows1 (n + 1) xs 
               | otherwise = rows1 (n + 1) xs
              
rowz l = rows1 0 l
            
split1 _ [] = []
split1 n l = (take n l) : (split1 n $ drop n l)

permuteLL a = rowlist $ map rowz $ transpose $ split1 b $ elems a
                 where b = fst $ snd $ bounds a

permute a | l == [] = error "singular matrix detected"
          | otherwise = tt
            where tt = l!!(head $ findIndices ( == (maximum n)) n)
                  n = map zz (dd a l)
                  zz a = foldr1 (+) (map abs (elems (diag a)))
                  dd a d | d == [] = [] 
                         | otherwise = mswap (head d) a : dd a (tail d)
                  l = permuteLL a
                 

rowlist :: Eq a => [[a]] -> [[a]]
rowlist [] = [[]]
rowlist (set:sets) = [x:xs | xs <- rowlist sets, x  <- set, not(x `elem` xs)]

diag :: Array (Int,Int) Double -- ^ A
     -> Array (Int,Int) Double -- ^ LU(A)

diag a = a'
    where a' = array bnds [ ((i,j), luij i j) | (i,j) <- range bnds ]
	  luij i j | i == j = a!(i,j)
	           | otherwise = 0
	  bnds = bounds a