diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,18 @@
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included
+in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+DEALINGS IN THE SOFTWARE.
+
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,4 @@
+import Distribution.Simple
+
+main = defaultMain
+
diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal
new file mode 100644
--- /dev/null
+++ b/bed-and-breakfast.cabal
@@ -0,0 +1,21 @@
+Name:           bed-and-breakfast
+Version:        0.1
+Synopsis:       Efficient Matrix operations in 100% Haskell.
+Description:    Efficient Matrix operations in 100% Haskell.
+                
+License:        MIT
+License-File:   LICENSE
+Author:         Julian Fleischer <julian.fleischer@fu-berlin.de>
+Maintainer:     Julian Fleischer <julian.fleischer@fu-berlin.de>
+Build-Type:     Simple
+Cabal-Version:  >= 1.8
+Category:       Data
+Stability:      stable
+
+Library
+    Exposed-Modules:    Numeric.Matrix
+    Build-Depends:      base >= 4.5 && < 5,
+                        array >= 0.4
+    Hs-Source-Dirs:     src
+
+
diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Matrix.hs
@@ -0,0 +1,513 @@
+{-# LANGUAGE Haskell2010
+    , TypeFamilies
+    , FlexibleContexts
+    , Trustworthy
+ #-}
+{-# OPTIONS -Wall -fno-warn-name-shadowing #-}
+
+module Numeric.Matrix (
+    Matrix,
+    MatrixElement (
+        matrix,
+        fromList,
+
+        unit,
+        zero,
+        diag,
+        empty,
+
+        at,
+        row,
+        col,
+
+        select,
+        toList,
+
+        dimensions,
+        numRows,
+        numCols,
+
+        isUnit,
+        isZero,
+        isDiagonal,
+        isEmpty,
+        isSquare,
+        
+        det,
+        rank,
+        transpose,
+        trace,
+
+        minus,
+        plus,
+        times,
+        inv,
+
+        map,
+        all,
+        any,
+
+        mapWithIndex,
+        allWithIndex,
+        anyWithIndex
+    )
+) where
+
+
+import Control.Applicative
+import Control.Monad
+import Control.Monad.ST
+
+import Data.Ratio
+import Data.Complex
+import qualified Data.List as L
+import Data.Array.IArray
+import Data.Array.MArray
+import Data.Array.Unboxed
+import Data.Array.ST
+import Data.STRef
+
+import Prelude hiding (any, all, read)
+import qualified Prelude as P
+
+
+data family Matrix e
+
+data instance Matrix Int    = IntMatrix    Int Int (Array Int (UArray Int Int))
+data instance Matrix Float  = FloatMatrix  Int Int (Array Int (UArray Int Float))
+data instance Matrix Double = DoubleMatrix Int Int (Array Int (UArray Int Double))
+
+data instance Matrix Integer = IntegerMatrix Int Int (Array Int (Array Int Integer))
+data instance Matrix (Ratio a) = RatioMatrix Int Int (Array Int (Array Int (Ratio a)))
+data instance Matrix (Complex a) = ComplexMatrix Int Int (Array Int (Array Int (Complex a)))
+
+instance (MatrixElement e, Show e) => Show (Matrix e) where
+    show = unlines . P.map showRow . toList
+      where
+        showRow = unwords . P.map ((' ':) . show)
+
+
+class Division e where
+    divide :: e -> e -> e
+
+instance Division Int    where divide = quot
+-- instance Division Int8   where divide = quot
+-- instance Division Int16  where divide = quot
+-- instance Division Int32  where divide = quot
+-- instance Division Int64  where divide = quot
+instance Division Integer where divide = quot
+instance Division Float  where divide = (/)
+instance Division Double where divide = (/)
+instance Integral a => Division (Ratio a) where divide = (/)
+instance RealFloat a => Division (Complex a) where divide = (/)
+
+
+
+class (Eq e, Num e) => MatrixElement e where
+
+    matrix :: (Int, Int) -> ((Int, Int) -> e) -> Matrix e
+    select :: ((Int, Int) -> Bool) -> Matrix e -> [e]
+    at :: Matrix e -> (Int, Int) -> e
+
+    row :: Int -> Matrix e -> [e]
+    col :: Int -> Matrix e -> [e]
+
+    dimensions :: Matrix e -> (Int, Int)
+    numRows :: Matrix e -> Int
+    numCols :: Matrix e -> Int
+
+    fromList :: [[e]] -> Matrix e
+    toList   :: Matrix e -> [[e]]
+
+    unit  :: Int -> Matrix e
+    zero  :: Int -> Matrix e
+    diag  :: [e] -> Matrix e
+    empty :: Matrix e
+
+    minus :: Matrix e -> Matrix e -> Matrix e
+    plus  :: Matrix e -> Matrix e -> Matrix e
+    times :: Matrix e -> Matrix e -> Matrix e
+    inv   :: Matrix e -> Maybe (Matrix e)
+
+--    adjugate  :: Matrix e -> Matrix e
+--    cofactors :: Matrix e -> Matrix e ; cofactors = undefined
+    det       :: Matrix e -> e
+    transpose :: Matrix e -> Matrix e
+    rank      :: Matrix e -> e
+    trace     :: Matrix e -> [e]
+
+    isUnit     :: Matrix e -> Bool
+    isDiagonal :: Matrix e -> Bool
+    isZero     :: Matrix e -> Bool
+    isEmpty    :: Matrix e -> Bool
+    isSquare   :: Matrix e -> Bool
+
+    select p m = [ at m (i,j) | i <- [1..numRows m]
+                              , j <- [1..numCols m]
+                              , p (i,j) ]
+
+    at m (i, j) = ((!! j) . (!! i) . toList) m
+
+    map :: MatrixElement f => (e -> f) -> Matrix e -> Matrix f
+    all :: (e -> Bool) -> Matrix e -> Bool
+    any :: (e -> Bool) -> Matrix e -> Bool
+
+    mapWithIndex :: MatrixElement f => ((Int, Int) -> e -> f) -> Matrix e -> Matrix f
+    allWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool
+    anyWithIndex :: ((Int, Int) -> e -> Bool) -> Matrix e -> Bool
+
+    unit n  = fromList [[ if i == j then 1 else 0 | j <- [1..n]] | i <- [1..n] ]
+    zero n  = matrix (n,n) (const 0)
+    empty   = fromList []
+    diag xs = matrix (n,n) (\(i,j) -> if i == j then xs !! (i-1) else 0)
+      where n = length xs
+    
+    row i m = ((!! (i-1)) . toList) m
+    col i m = (row i . transpose) m
+
+    numRows = fst . dimensions
+    numCols = snd . dimensions
+    dimensions m = case toList m of [] -> (0, 0)
+                                    (x:xs) -> (length xs + 1, length x)
+
+--    adjugate = transpose . cofactors
+    transpose = fromList . L.transpose . toList
+    trace = select (uncurry (==))
+    inv _ = Nothing
+
+    isZero = all (== 0)
+    isUnit m = isSquare m && P.all (== 1) (trace m)
+    isEmpty m = numRows m == 0 || numCols m == 0
+    isDiagonal = allWithIndex (uncurry $ \x y z -> if x /= y then z == 0 else True)
+    isSquare m = let (a, b) = dimensions m in a == b
+
+    map f = mapWithIndex (const f)
+    all f = allWithIndex (const f)
+    any f = anyWithIndex (const f)
+
+    mapWithIndex f m = matrix (dimensions m) (\x -> f x (m `at` x))
+    allWithIndex f m = P.all id [ f (i, j) (m `at` (i,j))
+                                | i <- [1..numRows m], j <- [1..numCols m]]
+    anyWithIndex f m = P.any id [ f (i, j) (m `at` (i,j))
+                                | i <- [1..numRows m], j <- [1..numCols m]]
+
+    a `plus` b
+        | dimensions a /= dimensions b = error "Matrix.plus: dimensions don't match."
+        | otherwise = matrix (dimensions a) (\x -> a `at` x + b `at` x)
+    a `minus` b
+        | dimensions a /= dimensions b = error "Matrix.minus: dimensions don't match."
+        | otherwise = matrix (dimensions a) (\x -> a `at` x - b `at` x)
+    a `times` b
+        | numRows a /= numCols b = error "Matrix.times: `numRows a' and `numCols b' don't match."
+        | otherwise = fromList [ [ row i a `dotProd` col j b
+                                 | j <- [1..numCols b] ]
+                               | i <- [1..numRows a] ]
+
+
+instance MatrixElement Int where
+    matrix   = _matrix IntMatrix
+    fromList = _fromList IntMatrix
+
+    at         (IntMatrix _ _ arr) = _at arr
+    dimensions (IntMatrix m n _) = (m, n)
+    row i      (IntMatrix _ _ arr) = _row i arr
+    col j      (IntMatrix _ _ arr) = _col j arr
+    toList     (IntMatrix _ _ arr) = _toList arr
+    inv = undefined -- IntMatrix $ runST (invSTU arr)
+    det        (IntMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)
+    rank = undefined -- runST (_rank thawsBoxed arr)
+
+instance MatrixElement Integer where
+    matrix   = _matrix IntegerMatrix
+    fromList = _fromList IntegerMatrix
+
+    at         (IntegerMatrix _ _ arr) = _at arr
+    dimensions (IntegerMatrix m n _) = (m, n)
+    row i      (IntegerMatrix _ _ arr) = _row i arr
+    col j      (IntegerMatrix _ _ arr) = _col j arr
+    toList     (IntegerMatrix _ _ arr) = _toList arr
+    inv = undefined -- IntMatrix $ runST (invSTU arr)
+    det        (IntegerMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)
+    rank = undefined -- runST (_rank thawsBoxed arr)
+
+instance MatrixElement Float where
+    matrix   = _matrix FloatMatrix
+    fromList = _fromList FloatMatrix
+
+    at         (FloatMatrix _ _ arr) = _at arr
+    dimensions (FloatMatrix m n _  ) = (m, n)
+    row i      (FloatMatrix _ _ arr) = _row i arr
+    col j      (FloatMatrix _ _ arr) = _col j arr
+    toList     (FloatMatrix _ _ arr) = _toList arr
+    inv        (FloatMatrix m n arr) = if m /= n then Nothing else
+                                        Just $ FloatMatrix m n $ runST (_inv unboxedST arr)
+    det        (FloatMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)
+    rank       (FloatMatrix _ _ arr) = runST (_rank thawsBoxed arr)
+
+instance MatrixElement Double where
+    matrix   = _matrix DoubleMatrix
+    fromList = _fromList DoubleMatrix
+
+    at         (DoubleMatrix _ _ arr) = _at arr
+    dimensions (DoubleMatrix m n _  ) = (m, n)
+    row i      (DoubleMatrix _ _ arr) = _row i arr
+    col j      (DoubleMatrix _ _ arr) = _col j arr
+    toList     (DoubleMatrix _ _ arr) = _toList arr
+    inv        (DoubleMatrix m n arr) = if m /= n then Nothing else
+                                         Just $ DoubleMatrix m n $ runST (_inv unboxedST arr)
+    det        (DoubleMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr)
+    rank       (DoubleMatrix _ _ arr) = runST (_rank thawsBoxed arr)
+
+instance Integral a => MatrixElement (Ratio a) where
+    matrix   = _matrix RatioMatrix
+    fromList = _fromList RatioMatrix
+
+    at         (RatioMatrix _ _ arr) = _at arr
+    dimensions (RatioMatrix m n _  ) = (m, n)
+    row i      (RatioMatrix _ _ arr) = _row i arr
+    col j      (RatioMatrix _ _ arr) = _col j arr
+    toList     (RatioMatrix _ _ arr) = _toList arr
+    inv        (RatioMatrix m n arr) = if m /= n then Nothing else
+                                        Just $ RatioMatrix m n $ runST (_inv boxedST arr)
+    det        (RatioMatrix m n arr) = if m /= n then 0 else  runST (_det thawsBoxed arr)
+    rank       (RatioMatrix _ _ arr) = runST (_rank thawsBoxed arr)
+
+instance RealFloat a => MatrixElement (Complex a) where
+    matrix   = _matrix ComplexMatrix
+    fromList = _fromList ComplexMatrix
+
+    at         (ComplexMatrix _ _ arr) = _at arr
+    dimensions (ComplexMatrix m n _  ) = (m, n)
+    row i      (ComplexMatrix _ _ arr) = _row i arr
+    col j      (ComplexMatrix _ _ arr) = _col j arr
+    toList     (ComplexMatrix _ _ arr) = _toList arr
+    inv        (ComplexMatrix m n arr) = if m /= n then Nothing else
+                                          Just $ ComplexMatrix m n $ runST (_inv boxedST arr)
+    det        (ComplexMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr)
+    rank       (ComplexMatrix _ _ arr) = runST (_rank thawsBoxed arr)
+
+
+_at :: (IArray a (u Int e), IArray u e)
+    => a Int (u Int e) -> (Int, Int) -> e
+_at arr (i,j) = arr ! i ! j
+
+_row, _col :: (IArray a (u Int e), IArray u e) => Int -> a Int (u Int e) -> [e]
+_row i arr = let row = arr ! i in [ row ! j | j <- [1..(snd (bounds arr))] ]
+_col j arr = [ arr ! i ! j | i <- [1..(snd (bounds arr))] ]
+
+_matrix :: (IArray a (u Int e), IArray u e)
+        => (Int -> Int -> (a Int (u Int e)) -> matrix e)
+        -> (Int, Int)
+        -> ((Int, Int) -> e)
+        -> matrix e
+_matrix c (numRows, numCols) generator =
+    c numRows numCols
+      $ array (1, numRows)
+      $ [ (i, array (1, numCols) [(j, generator (i, j))
+                                 | j <- [1..numCols]])
+          | i <- [1..numRows] ]
+
+_toList :: (IArray a e) => Array Int (a Int e) -> [[e]]
+_toList = P.map elems . elems
+
+_fromList :: (IArray a (u Int e), IArray u e)
+          => (Int -> Int -> a Int (u Int e) -> matrix e) -> [[e]] -> matrix e
+_fromList c xs =
+    let lengths = P.map length xs
+        numCols = foldl1 min lengths
+        numRows = length lengths
+        
+    in  c numRows numCols
+          $ array (1, numRows)
+          $ zip [1..numRows]
+          $ P.map (array (1, numCols) . zip [1..numCols]) xs
+
+dotProd :: Num a => [a] -> [a] -> a
+dotProd x = L.foldl' (+) 0 . zipWith (*) x
+
+thawsBoxed :: (IArray a e, MArray (STArray s) e (ST s))
+           => Array Int (a Int e)
+           -> ST s [STArray s Int e]
+thawsBoxed = mapM thaw . elems
+
+thawsUnboxed :: (IArray a e, MArray (STUArray s) e (ST s))
+             => Array Int (a Int e)
+             -> ST s [STUArray s Int e]
+thawsUnboxed = mapM thaw . elems
+
+arrays :: [(u s) Int e]
+       -> ST s ((STArray s) Int ((u s) Int e))
+arrays list = newListArray (1, length list) list
+
+augment :: (IArray a e, MArray (u s) e (ST s), Num e)
+        => ((Int, Int) -> [e] -> ST s ((u s) Int e))
+        -> Array Int (a Int e)
+        -> ST s (STArray s Int (u s Int e))
+augment _ arr = do
+    let (_, n) = bounds arr
+        row (a,i) = newListArray (1, 2*n)
+                                 [ if j > n then (if j == i + n then 1 else 0)
+                                            else a ! j
+                                 | j <- [1..(2*n)] ]
+    
+    mapM row (zip (elems arr) [1..]) >>= newListArray (1, n)
+
+boxedST :: MArray (STArray s) e (ST s)
+        => (Int, Int) -> [e] -> ST s ((STArray s) Int e)
+boxedST = newListArray
+
+unboxedST :: MArray (STUArray s) e (ST s)
+          => (Int, Int) -> [e] -> ST s ((STUArray s) Int e)
+unboxedST = newListArray
+
+
+tee :: Monad m => (b -> m a) -> b -> m b
+tee f x = f x >> return x
+
+read :: (MArray a1 b m, MArray a (a1 Int b) m) =>
+                       a Int (a1 Int b) -> Int -> Int -> m b
+read a i j = readArray a i >>= flip readArray j
+
+
+_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e)
+     => ((Int, Int) -> [e] -> ST s ((u s) Int e))
+     -> Array Int (a Int e)
+     -> ST s (Array Int (a Int e))
+_inv mkArrayST mat = do
+    let m = snd $ bounds mat
+        n = 2*m
+
+    a <- augment mkArrayST mat
+
+    flip mapM_ [1..m] $ \k -> do
+        flip mapM_ [(k+1)..m] $ \i -> do
+            a_i <- readArray a i
+            a_k <- readArray a k
+            flip mapM_ [(k+1)..n] $ \j -> do
+                a_ij <- readArray a_i j
+                a_kj <- readArray a_k j
+                a_ik <- readArray a_i k
+                a_kk <- readArray a_k k
+                writeArray a_i j (a_ij - a_kj * (a_ik / a_kk))
+            writeArray a_i k 0
+
+    flip mapM_ [ m - k | k <- [1..(m-1)] ] $ \i -> do
+        r1 <- readArray a i
+        r2 <- readArray a (i+1)
+
+        p <- readArray r2 (i+1) >>= return . (1 /)
+
+        flip mapM_ [(i+1)..2*m] $ \j -> do
+            c1 <- readArray r1 j
+            c2 <- readArray r2 j
+            writeArray r1 j (c2 * p + c1)
+
+    result <- flip mapM [1..m] $ \i -> do
+        r <- readArray a i
+        p <- readArray r i
+        
+        mapM (\j -> (/ p) <$> readArray r j) [(m+1)..(2*m)]
+            >>= return . listArray (1, m)
+    
+    return $ listArray (1, m) result
+
+
+_rank :: (IArray a e, MArray (u s) e (ST s),
+           Fractional e, Eq e)
+      => (Array Int (a Int e) -> ST s [(u s) Int e])
+      -> Array Int (a Int e)
+      -> ST s e
+_rank thaws mat = do
+    let m = snd $ bounds mat
+        n = snd $ bounds (mat ! 1)
+
+    a <- thaws mat >>= arrays
+
+    trace <- flip mapM [1..m] $ \k -> do
+        flip mapM_ [(k+1)..m] $ \i -> do
+            a_i <- readArray a i
+            a_k <- readArray a k
+            flip mapM_ [(k+1)..n] $ \j -> do
+                a_ij <- readArray a_i j
+                a_kj <- readArray a_k j
+                a_ik <- readArray a_i k
+                a_kk <- readArray a_k k
+                writeArray a_i j (a_ij - a_kj * (a_ik / a_kk))
+            writeArray a_i k 0
+        read a k k
+
+    return $ fromIntegral $ length $ filter (/= 0) trace
+
+
+_det :: (IArray a e, MArray (u s) e (ST s),
+         Num e, Eq e, Division e)
+     => (Array Int (a Int e) -> ST s [(u s) Int e])
+     -> Array Int (a Int e) -> ST s e
+_det thaws mat = do
+
+    let size = snd $ bounds mat
+
+    a <- thaws mat >>= arrays
+
+    signR  <- newSTRef 1
+    pivotR <- newSTRef 1
+
+    flip mapM_ [1..size] $ \k -> do
+
+        prev  <- readSTRef pivotR
+        pivot <- read a k k >>= tee (writeSTRef pivotR)
+
+        when (pivot == 0) $ do
+            s <- flip mapM [(k+1)..size] $ \r -> do
+                a_rk <- read a r k
+                if a_rk == 0 then return 0 else return r
+            let sf = filter (>0) s
+
+            when (not $ null sf) $ do
+                let sw = head sf
+
+                row <- readArray a sw
+                readArray a k >>= writeArray a sw
+                writeArray a k row
+
+                read a k k >>= writeSTRef pivotR
+                readSTRef signR >>= writeSTRef signR . negate
+
+        pivot' <- readSTRef pivotR
+        flip mapM [(k+1)..size] $ \i -> do
+            a_i <- readArray a i
+            flip mapM [(k+1)..size] $ \j -> do
+                a_ij <- readArray a_i j
+                a_ik <- readArray a_i k
+                a_kj <- read a k j
+                writeArray a_i j ((pivot' * a_ij - a_ik * a_kj) `divide` prev)
+
+    liftM2 (*) (readSTRef pivotR) (readSTRef signR)
+
+
+mat n = fromList [ [ j | j <- take n [i,(i^i)..] ] | i <- take n [1..] ]
+
+
+{-
+-- | The 'findIndex' function takes a predicate and a list and returns
+-- the index of the first element in the list satisfying the predicate,
+-- or 'Nothing' if there is no such element.
+findIndex       :: (a -> Bool) -> [a] -> Maybe Int
+findIndex p     = listToMaybe . findIndices p
+
+-- | The 'findIndices' function extends 'findIndex', by returning the
+-- indices of all elements satisfying the predicate, in ascending order.
+findIndices      :: (a -> Bool) -> [a] -> [Int]
+
+#if defined(USE_REPORT_PRELUDE) || !defined(__GLASGOW_HASKELL__)
+findIndices p xs = [ i | (x,i) <- zip xs [0..], p x]
+#else
+-- Efficient definition
+findIndices p ls = loop 0# ls
+                 where
+                   loop _ [] = []
+                   loop n (x:xs) | p x       = I# n : loop (n +# 1#) xs
+                                 | otherwise = loop (n +# 1#) xs
+#endif  /* USE_REPORT_PRELUDE */
+-}
+