diff --git a/Data/._.DS_Store b/Data/._.DS_Store
deleted file mode 100644
Binary files a/Data/._.DS_Store and /dev/null differ
diff --git a/Data/Eigen/Internal.hi b/Data/Eigen/Internal.hi
new file mode 100644
Binary files /dev/null and b/Data/Eigen/Internal.hi differ
diff --git a/Data/Eigen/Internal.hs b/Data/Eigen/Internal.hs
--- a/Data/Eigen/Internal.hs
+++ b/Data/Eigen/Internal.hs
@@ -5,8 +5,7 @@
 import Foreign.C.Types
 import Foreign.C.String
 import Control.Monad
-
-data C_MatrixXd
+import System.IO.Unsafe
 
 class Cast a b where
     cast :: a -> b
@@ -23,34 +22,28 @@
 instance Cast Int CInt where
     cast = fromIntegral
 
+performIO :: IO a -> a
+performIO = unsafeDupablePerformIO
+
 foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()
 
 call :: IO CString -> IO ()
 call func = func >>= \c_str -> when (c_str /= nullPtr) $
-    peekCString c_str >>= \str -> c_freeString c_str >> error str
+    peekCString c_str >>= \str -> c_freeString c_str >> fail str
 
 foreign import ccall "eigen-proxy.h eigen_initParallel" c_initParallel :: IO ()
 foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO ()
 
-foreign import ccall "eigen-proxy.h eigen_create"      c_create :: CInt -> CInt -> IO (Ptr C_MatrixXd)
-foreign import ccall "eigen-proxy.h &eigen_destroy"    c_destroy :: FunPtr (Ptr C_MatrixXd -> IO ())
-foreign import ccall "eigen-proxy.h eigen_clone"       c_clone :: Ptr C_MatrixXd -> IO (Ptr C_MatrixXd)
-foreign import ccall "eigen-proxy.h eigen_get"         c_get :: Ptr CDouble -> Ptr C_MatrixXd -> CInt -> CInt -> IO CString
-foreign import ccall "eigen-proxy.h eigen_set"         c_set :: Ptr C_MatrixXd -> CInt -> CInt -> CDouble -> IO CString
-foreign import ccall "eigen-proxy.h eigen_data"        c_data :: Ptr C_MatrixXd -> IO (Ptr CDouble)
-foreign import ccall "eigen-proxy.h eigen_rows"        c_rows :: Ptr C_MatrixXd -> IO CInt
-foreign import ccall "eigen-proxy.h eigen_cols"        c_cols :: Ptr C_MatrixXd -> IO CInt
-foreign import ccall "eigen-proxy.h eigen_copy"        c_copy :: Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO ()
-foreign import ccall "eigen-proxy.h eigen_resize"      c_resize :: Ptr C_MatrixXd -> CInt -> CInt -> IO ()
-foreign import ccall "eigen-proxy.h eigen_add"         c_add :: Ptr C_MatrixXd -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_sub"         c_sub :: Ptr C_MatrixXd -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_mul"         c_mul :: Ptr C_MatrixXd -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_transpose"   c_transpose :: Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_inverse"     c_inverse :: Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_adjoint"     c_adjoint :: Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_normalize"   c_normalize :: Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_norm"        c_norm :: Ptr C_MatrixXd -> IO CDouble
-foreign import ccall "eigen-proxy.h eigen_squaredNorm" c_squaredNorm :: Ptr C_MatrixXd -> IO CDouble
-foreign import ccall "eigen-proxy.h eigen_blueNorm"    c_blueNorm :: Ptr C_MatrixXd -> IO CDouble
-foreign import ccall "eigen-proxy.h eigen_hypotNorm"   c_hypotNorm :: Ptr C_MatrixXd -> IO CDouble
-foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: Ptr C_MatrixXd -> IO CDouble
+foreign import ccall "eigen-proxy.h eigen_add"         c_add         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_sub"         c_sub         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_mul"         c_mul         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_transpose"   c_transpose   :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_inverse"     c_inverse     :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_adjoint"     c_adjoint     :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_conjugate"   c_conjugate   :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_normalize"   c_normalize   :: Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_norm"        c_norm        :: Ptr CDouble -> CInt -> CInt -> IO CDouble
+foreign import ccall "eigen-proxy.h eigen_squaredNorm" c_squaredNorm :: Ptr CDouble -> CInt -> CInt -> IO CDouble
+foreign import ccall "eigen-proxy.h eigen_blueNorm"    c_blueNorm    :: Ptr CDouble -> CInt -> CInt -> IO CDouble
+foreign import ccall "eigen-proxy.h eigen_hypotNorm"   c_hypotNorm   :: Ptr CDouble -> CInt -> CInt -> IO CDouble
+foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: Ptr CDouble -> CInt -> CInt -> IO CDouble
diff --git a/Data/Eigen/Internal.o b/Data/Eigen/Internal.o
new file mode 100644
Binary files /dev/null and b/Data/Eigen/Internal.o differ
diff --git a/Data/Eigen/LA.hi b/Data/Eigen/LA.hi
new file mode 100644
Binary files /dev/null and b/Data/Eigen/LA.hi differ
diff --git a/Data/Eigen/LA.hs b/Data/Eigen/LA.hs
--- a/Data/Eigen/LA.hs
+++ b/Data/Eigen/LA.hs
@@ -66,13 +66,14 @@
 import Foreign.C.String
 import Foreign.Storable
 import Foreign.Marshal.Alloc
+import Control.Applicative
 import Data.Eigen.Matrix
-import qualified Data.Eigen.Matrix.Mutable as MM
 import Data.Eigen.Internal
-import System.IO.Unsafe
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Storable.Mutable as VSM
 
-foreign import ccall "eigen-proxy.h eigen_solve" c_solve :: CInt -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO CString
-foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO CString
+foreign import ccall "eigen-proxy.h eigen_solve" c_solve :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
+foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString
 
 
 
@@ -113,23 +114,35 @@
 
 -- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@
 solve :: Decomposition -> Matrix -> Matrix -> Matrix
-solve d a b = (`modify` empty) $ \x ->
-    with a $ \pa ->
-    with b $ \pb ->
-    MM.with x $ \px ->
-        call $ c_solve (fromIntegral $ fromEnum d) px pa pb
+solve d a b = performIO $ do
+    let
+        cols = 1
+        rows = m_cols a
+    vals <- VSM.new (rows * cols)
+    VSM.unsafeWith vals $ \px ->
+        VS.unsafeWith (m_vals a) $ \pa ->
+            VS.unsafeWith (m_vals b) $ \pb ->
+                call $ c_solve (cast $ fromEnum d)
+                    px (cast rows) (cast cols)
+                    pa (cast $ m_rows a) (cast $ m_cols a)
+                    pb (cast $ m_rows b) (cast $ m_cols b)
+    Matrix rows cols <$> VS.unsafeFreeze vals
 
+
 -- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm
 relativeError :: Matrix -> Matrix -> Matrix -> Double
-relativeError x a b = unsafePerformIO $
-    with x $ \px ->
-    with a $ \pa ->
-    with b $ \pb ->
-    alloca $ \pr -> do
-        call $ c_relativeError pr px pa pb
-        fmap cast $ peek pr
+relativeError x a b = performIO $ do
+    VS.unsafeWith (m_vals x) $ \px ->
+        VS.unsafeWith (m_vals a) $ \pa ->
+            VS.unsafeWith (m_vals b) $ \pb ->
+                alloca $ \pe -> do
+                    call $ c_relativeError pe
+                        px (cast $ m_rows x) (cast $ m_cols x)
+                        pa (cast $ m_rows a) (cast $ m_cols a)
+                        pb (cast $ m_rows b) (cast $ m_cols b)
+                    cast <$> peek pe
 
-{- | 
+{- |
 [(coeffs, error) = linearRegression points] computes multiple linear regression @y = a1 x1 + a2 x2 + ... + an xn + b@ using 'ColPivHouseholderQR' decomposition
 
 * point format is @[y, x1..xn]@
@@ -153,7 +166,7 @@
  @
  ([-2.3466569233817127,-0.2534897541434826,-0.1749653335680988],1.8905965120153139e-3)
  @
- 
+
 -}
 linearRegression :: [[Double]] -> ([Double], Double)
 linearRegression points = (coeffs, e) where
diff --git a/Data/Eigen/LA.o b/Data/Eigen/LA.o
new file mode 100644
Binary files /dev/null and b/Data/Eigen/LA.o differ
diff --git a/Data/Eigen/Matrix.hi b/Data/Eigen/Matrix.hi
new file mode 100644
Binary files /dev/null and b/Data/Eigen/Matrix.hi differ
diff --git a/Data/Eigen/Matrix.hs b/Data/Eigen/Matrix.hs
--- a/Data/Eigen/Matrix.hs
+++ b/Data/Eigen/Matrix.hs
@@ -1,10 +1,11 @@
-{-# LANGUAGE RecordWildCards, MultiParamTypeClasses #-}
+{-# LANGUAGE RecordWildCards, MultiParamTypeClasses, Rank2Types #-}
 module Data.Eigen.Matrix (
     -- * Matrix type
     Matrix(..),
     -- * Matrix conversions
     fromList,
     toList,
+    generate,
     -- * Standard matrices and special cases
     empty,
     zero,
@@ -28,32 +29,37 @@
     norm,
     squaredNorm,
     determinant,
+    -- * Matrix operations
+    add,
+    sub,
+    mul,
     -- * Matrix transformations
     inverse,
     adjoint,
+    conjugate,
     transpose,
     normalize,
-    -- * Mutable operations
-    freeze,
-    thaw,
     modify,
-    with
+    -- * Mutable matrices
+    thaw,
+    freeze,
+    unsafeThaw,
+    unsafeFreeze
 ) where
 
 import Data.List (intercalate)
-import Text.Printf
+import Data.Tuple
 import Foreign.Ptr
-import Foreign.Storable
-import Foreign.ForeignPtr
 import Foreign.C.Types
-import Foreign.Marshal.Array
+import Foreign.C.String
 import Control.Monad
 import Control.Monad.ST
-import System.IO.Unsafe
+import Control.Monad.Primitive
+import Control.Applicative hiding (empty)
 import qualified Data.Vector.Storable as VS
 import qualified Data.Vector.Storable.Mutable as VSM
-import qualified Data.Eigen.Matrix.Mutable as MM
 import Data.Eigen.Internal
+import Data.Eigen.Matrix.Mutable
 
 -- | constant Matrix class to be used in pure computations, uses the same column major memory layout as Eigen MatrixXd
 data Matrix = Matrix {
@@ -64,14 +70,14 @@
 
 -- | pretty prints the matrix
 instance Show Matrix where
-    show m@Matrix{..} = printf "Matrix %dx%d\n%s\n" m_rows m_cols $
-        intercalate "\n" $ map (intercalate "\t" . map show) $ toList m
+    show m@Matrix{..} = concat [
+        "Matrix ", show m_rows, "x", show m_cols, "\n", intercalate "\n" $ map (intercalate "\t" . map show) $ toList m, "\n"]
 
 -- | only the following functions are defined for Num instance: (*), (+), (-)
 instance Num Matrix where
-    (*) = binop MM.mul
-    (+) = binop MM.add
-    (-) = binop MM.sub
+    (*) = mul
+    (+) = add
+    (-) = sub
     fromInteger = undefined
     signum = undefined
     abs = undefined
@@ -94,11 +100,11 @@
 
 -- | square matrix with 1 on main diagonal and 0 elsewhere
 identity :: Int -> Matrix
-identity size = Matrix size size $ runST $ do
+identity size = Matrix size size $ VS.create $ do
     vm <- VSM.replicate (size * size) 0
     forM_ [0..pred size] $ \n ->
         VSM.write vm (n * size + n) 1
-    VS.unsafeFreeze vm
+    return vm
 
 -- | number of rows for the matrix
 rows :: Matrix -> Int
@@ -122,8 +128,9 @@
 
 -- | extract rectangular block from matrix defined by startRow startCol blockRows blockCols
 block :: Int -> Int -> Int -> Int -> Matrix -> Matrix
-block startRow startCol blockRows blockCols m = fromList $
-    [[coeff row col m | col <- take blockCols [startCol..]] | row <- take blockRows [startRow..]]
+block startRow startCol blockRows blockCols m =
+    generate blockRows blockCols $ \row col ->
+        coeff (startRow + row) (startCol + col) m
 
 -- | the maximum of all coefficients of matrix
 maxCoeff :: Matrix -> Double
@@ -149,90 +156,130 @@
 rightCols :: Int -> Matrix -> Matrix
 rightCols cols m@Matrix{..} = block 0 (m_cols - cols) m_rows cols m
 
-binop :: (MM.MMatrix -> MM.MMatrix -> MM.MMatrix -> IO a) -> Matrix -> Matrix -> Matrix
-binop f lhs rhs = unsafePerformIO $ do
-    ret <- MM.new 0 0
-    lhs <- thaw lhs
-    rhs <- thaw rhs
-    _ <- f ret lhs rhs
-    freeze ret
-
 -- | construct matrix from a list of rows, column count is detected as maximum row length
 fromList :: [[Double]] -> Matrix
 fromList list = Matrix rows cols vals where
     rows = length list
     cols = maximum $ map length list
-    vals = runST $ do
+    vals = VS.create $ do
         vm <- VSM.replicate (rows * cols) 0
-        zipWithM_ (\row vals ->
-            zipWithM_ (\col val ->
-                VSM.write vm (col * rows + row) (cast val)) [0..] vals) [0..] list
-        VS.unsafeFreeze vm
+        forM_ (zip [0..] list) $ \(row, vals) ->
+            forM_ (zip [0..] vals) $ \(col, val) ->
+                VSM.write vm (col * rows + row) (cast val)
+        return vm
 
 -- | converts matrix to a list of its rows
 toList :: Matrix -> [[Double]]
 toList Matrix{..} = [[cast $ m_vals VS.! (col * m_rows + row) | col <- [0..pred m_cols]] | row <- [0..pred m_rows]]
 
+-- | craete matrix using generator function f :: row -> col -> val
+generate :: Int -> Int -> (Int -> Int -> Double) -> Matrix
+generate rows cols f = Matrix rows cols $ VS.create $ do
+    vals <- VSM.new (rows * cols)
+    forM_ [0..pred rows] $ \row ->
+        forM_ [0..pred cols] $ \col ->
+            VSM.write vals (col * rows + row) (cast $ f row col)
+    return vals
+
+
 -- | for vectors, the l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the square root of the sum of the square of all the matrix entries. For vectors, this is also equals to the square root of the dot product of this with itself.
 norm :: Matrix -> Double
-norm m = unsafePerformIO $ thaw m >>= MM.norm
+norm = _unop c_norm
 
 -- | for vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.
 squaredNorm :: Matrix -> Double
-squaredNorm m = unsafePerformIO $ thaw m >>= MM.squaredNorm
+squaredNorm = _unop c_squaredNorm
 
 -- | the determinant of the matrix
 determinant :: Matrix -> Double
-determinant m = unsafePerformIO $ thaw m >>= MM.determinant
+determinant m@Matrix{..}
+    | m_cols == m_rows = _unop c_determinant m
+    | otherwise = error "you tried calling determinant on non-square matrix"
 
+-- | return a - b
+add :: Matrix -> Matrix -> Matrix
+add = _binop c_add
+
+-- | return a + b
+sub :: Matrix -> Matrix -> Matrix
+sub = _binop c_sub
+
+-- | return a * b
+mul :: Matrix -> Matrix -> Matrix
+mul = _binop c_mul
+
+
 {- | inverse of the matrix
 
 For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses class 'PartialPivLU'
 -}
 inverse :: Matrix -> Matrix
-inverse = modify MM.inverse
+inverse m@Matrix{..}
+    | m_rows == m_cols = _modify id c_inverse m
+    | otherwise = error "you tried calling inverse on non-square matrix"
 
 -- | adjoint of the matrix
 adjoint :: Matrix -> Matrix
-adjoint = modify MM.adjoint
+adjoint = _modify swap c_adjoint
 
 -- | transpose of the matrix
 transpose :: Matrix -> Matrix
-transpose = modify MM.transpose
+transpose = _modify swap c_transpose
 
+-- | conjugate of the matrix
+conjugate :: Matrix -> Matrix
+conjugate = _modify id c_conjugate
+
 -- | nomalize the matrix by deviding it on its 'norm'
 normalize :: Matrix -> Matrix
-normalize = modify MM.normalize
+normalize Matrix{..} = performIO $ do
+    vals <- VS.thaw m_vals
+    VSM.unsafeWith vals $ \p ->
+        call $ c_normalize p (cast m_rows) (cast m_cols)
+    Matrix m_rows m_cols <$> VS.unsafeFreeze vals
 
--- | create a snapshot of mutable matrix
-freeze :: MM.MMatrix -> IO Matrix
-freeze mm = MM.with mm $ \pm -> do
-    rows <- fmap cast $ c_rows pm
-    cols <- fmap cast $ c_cols pm
-    let len = rows * cols
-    src <- c_data pm
-    fp <- mallocForeignPtrArray len
-    withForeignPtr fp $ \dst -> copyArray dst src len
-    return Matrix {
-        m_rows = rows,
-        m_cols = cols,
-        m_vals = VS.unsafeFromForeignPtr fp 0 len
-    }
+-- | Apply a destructive operation to a matrix. The operation will be performed in place if it is safe to do so and will modify a copy of the matrix otherwise.
+modify :: (forall s. MMatrix s -> ST s ()) -> Matrix -> Matrix
+modify f m@Matrix{..} = m { m_vals = VS.modify f' m_vals } where
+    f' vals = f (MMatrix m_rows m_cols vals)
 
--- | create mutable copy of the matrix
-thaw :: Matrix -> IO MM.MMatrix
-thaw Matrix{..} = withForeignPtr fp $ \src -> do
-    let len = m_rows * m_cols
-    pm <- c_create (cast m_rows) (cast m_cols)
-    dst <- c_data pm
-    copyArray dst (plusPtr src $ off * sizeOf (undefined :: CDouble)) len
-    fmap MM.MMatrix $ newForeignPtr c_destroy pm
-    where (fp, off, _) = VS.unsafeToForeignPtr m_vals
+-- | Yield an immutable copy of the mutable matrix
+freeze :: PrimMonad m => MMatrix (PrimState m) -> m Matrix
+freeze MMatrix{..} = VS.freeze mm_vals >>= \vals -> return $ Matrix mm_rows mm_cols vals
 
--- | apply mutable operation to the mutable copy of the matrix and snapshot of this copy
-modify :: (MM.MMatrix -> IO ()) -> Matrix -> Matrix
-modify f m = unsafePerformIO $ thaw m >>= \mm -> f mm >> freeze mm
+-- | Yield a mutable copy of the immutable matrix
+thaw :: PrimMonad m => Matrix -> m (MMatrix (PrimState m))
+thaw Matrix{..} = VS.thaw m_vals >>= \vals -> return $ MMatrix m_rows m_cols vals
 
--- | apply foreign operation to the mutable copy of the matrix and operation result
-with :: Matrix -> (Ptr C_MatrixXd -> IO a) -> IO a
-with m f = thaw m >>= \mm -> MM.with mm f
+-- | Unsafe convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.
+unsafeFreeze :: PrimMonad m => MMatrix (PrimState m) -> m Matrix
+unsafeFreeze MMatrix{..} = VS.unsafeFreeze mm_vals >>= \vals -> return $ Matrix mm_rows mm_cols vals
+
+-- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.
+unsafeThaw :: PrimMonad m => Matrix -> m (MMatrix (PrimState m))
+unsafeThaw Matrix{..} = VS.unsafeThaw m_vals >>= \vals -> return $ MMatrix m_rows m_cols vals
+
+_unop :: (Ptr CDouble -> CInt -> CInt -> IO CDouble) -> Matrix -> Double
+_unop f Matrix{..} = performIO $ VS.unsafeWith m_vals $ \p ->
+    cast <$> f p (cast m_rows) (cast m_cols)
+
+_binop :: (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix -> Matrix
+_binop f m1 m2 = performIO $ do
+    vals <- VS.thaw (m_vals m1)
+    VSM.unsafeWith vals $ \lhs ->
+        VS.unsafeWith (m_vals m2) $ \rhs ->
+            call $ f
+                lhs (cast $ m_rows m1) (cast $ m_cols m1)
+                rhs (cast $ m_rows m2) (cast $ m_cols m2)
+    Matrix (m_rows m1) (m_cols m1) <$> VS.unsafeFreeze vals
+
+_modify :: ((Int,Int) -> (Int,Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix
+_modify f g Matrix{..} = performIO $ do
+    let (rows, cols) = f (m_rows, m_cols)
+    vals <- VSM.new (rows * cols)
+    VS.unsafeWith m_vals $ \src ->
+        VSM.unsafeWith vals $ \dst ->
+            call $ g
+                dst (cast rows) (cast cols)
+                src (cast m_rows) (cast m_cols)
+    Matrix rows cols <$> VS.unsafeFreeze vals
diff --git a/Data/Eigen/Matrix.o b/Data/Eigen/Matrix.o
new file mode 100644
Binary files /dev/null and b/Data/Eigen/Matrix.o differ
diff --git a/Data/Eigen/Matrix/Mutable.hi b/Data/Eigen/Matrix/Mutable.hi
new file mode 100644
Binary files /dev/null and b/Data/Eigen/Matrix/Mutable.hi differ
diff --git a/Data/Eigen/Matrix/Mutable.hs b/Data/Eigen/Matrix/Mutable.hs
--- a/Data/Eigen/Matrix/Mutable.hs
+++ b/Data/Eigen/Matrix/Mutable.hs
@@ -1,133 +1,46 @@
-{-# LANGUAGE RecordWildCards, MultiParamTypeClasses #-}
-module Data.Eigen.Matrix.Mutable (
-    MMatrix(..),
-    new,
-    rows,
-    cols,
-    copy,
-    clone,
-    resize,
-    get,
-    set,
-    setRow,
-    setCol,
-    add,
-    sub,
-    mul,
-    inverse,
-    adjoint,
-    transpose,
-    normalize,
-    norm,
-    blueNorm,
-    hypotNorm,
-    squaredNorm,
-    determinant,
-    with
-) where
+{-# LANGUAGE RecordWildCards #-}
+module Data.Eigen.Matrix.Mutable where
 
-import Foreign.Ptr
-import Foreign.ForeignPtr
-import Foreign.C.String
-import Foreign.Storable
-import Foreign.Marshal.Alloc
-import Control.Monad
+import Control.Monad.Primitive
+import Foreign.C.Types
+import qualified Data.Vector.Storable.Mutable as VSM
 import Data.Eigen.Internal
 
-data MMatrix = MMatrix { mm_fp :: ForeignPtr C_MatrixXd };
+-- | Mutable matrix. You can modify elements
+data MMatrix s = MMatrix {
+    mm_rows :: Int,
+    mm_cols :: Int,
+    mm_vals :: VSM.MVector s CDouble
+}
 
-new :: Int -> Int -> IO MMatrix
+type IOMatrix = MMatrix RealWorld
+type STMatrix s = MMatrix s
+
+-- | Creates a mutable matrix of the given dimension. Elements are initialized with 0.
+new :: PrimMonad m => Int -> Int -> m (MMatrix (PrimState m))
 new rows cols = do
-    pm <- c_create (cast rows) (cast cols)
-    fp <- newForeignPtr c_destroy pm
-    return $ MMatrix fp
+    vals <- VSM.replicate (rows * cols) 0
+    return $ MMatrix rows cols vals
 
-rows :: MMatrix -> IO Int
-rows mm = fmap cast $ with mm c_rows
+-- | Set all elements of the matrix to the given value
+set :: PrimMonad m => (MMatrix (PrimState m)) -> Double -> m ()
+set MMatrix{..} val = VSM.set mm_vals (cast val)
 
-cols :: MMatrix -> IO Int
-cols mm = fmap cast $ with mm c_cols
+-- | Yield the element at the given position.
+read :: PrimMonad m => MMatrix (PrimState m) -> Int -> Int -> m Double
+read MMatrix{..} row col
+    | row >= 0 && row < mm_rows && 
+      col >= 0 && col < mm_cols = do
+        val <- VSM.read mm_vals (col * mm_rows + row)
+        return $ cast val
+    | otherwise = fail "MMatrix.read: wrong index"
 
-copy :: MMatrix -> MMatrix -> IO ()
-copy dst src =
-    with dst $ \pdst ->
-    with src $ \psrc ->
-        c_copy pdst psrc
+-- | Replace the element at the given position.
+write :: PrimMonad m => MMatrix (PrimState m) -> Int -> Int -> Double -> m ()
+write MMatrix{..} row col val
+    | row >= 0 && row < mm_rows && 
+      col >= 0 && col < mm_cols = 
+        VSM.write mm_vals (col * mm_rows + row) (cast val)
+    | otherwise = fail "MMatrix.write: wrong index"
 
-clone :: MMatrix -> IO MMatrix
-clone mm = with mm $ \src -> do
-    dst <- c_clone src
-    fp <- newForeignPtr c_destroy dst
-    return $ MMatrix fp
 
-resize :: MMatrix -> Int -> Int -> IO ()
-resize mm rows cols = with mm $ \pm ->
-    c_resize pm (cast rows) (cast cols)
-
-get :: MMatrix -> Int -> Int -> IO Double
-get mm row col =
-    with mm $ \pm ->
-    alloca $ \pv -> do
-        call $ c_get pv pm (cast row) (cast col)
-        fmap cast $ peek pv
-
-set :: MMatrix -> Int -> Int -> Double -> IO ()
-set mm row col val = with mm $ \pm -> call $ c_set pm (cast row) (cast col) (cast val)
-
-setRow :: MMatrix -> Int -> [Double] -> IO ()
-setRow mm row vals = zipWithM_ (\col val -> set mm row col val) [0..] vals
-
-setCol :: MMatrix -> Int -> [Double] -> IO ()
-setCol mm col vals = zipWithM_ (\row val -> set mm row col val) [0..] vals
-
-binop :: (Ptr C_MatrixXd -> Ptr C_MatrixXd -> Ptr C_MatrixXd -> IO CString) -> MMatrix -> MMatrix -> MMatrix -> IO ()
-binop f ret lhs rhs  =
-    with ret $ \pret ->
-    with lhs $ \plhs ->
-    with rhs $ \prhs ->
-        call $ f pret plhs prhs
-
-add :: MMatrix -> MMatrix -> MMatrix -> IO ()
-add = binop c_add
-
-sub :: MMatrix -> MMatrix -> MMatrix -> IO ()
-sub = binop c_sub
-
-mul :: MMatrix -> MMatrix -> MMatrix -> IO ()
-mul = binop c_mul
-
-inpop :: (Ptr C_MatrixXd -> IO CString) -> MMatrix -> IO ()
-inpop f mm = with mm $ \pm -> call $ f pm
-
-inverse :: MMatrix -> IO ()
-inverse = inpop c_inverse
-
-adjoint :: MMatrix -> IO ()
-adjoint = inpop c_adjoint
-
-transpose :: MMatrix -> IO ()
-transpose = inpop c_transpose
-
-normalize :: MMatrix -> IO ()
-normalize = inpop c_normalize
-
-norm :: MMatrix -> IO Double
-norm mm = fmap cast $ with mm c_norm
-
-squaredNorm :: MMatrix -> IO Double
-squaredNorm mm = fmap cast $ with mm c_squaredNorm
-
---stableNorm :: MMatrix -> IO Double
---stableNorm mm = fmap cast $ with mm c_stableNorm
-
-blueNorm :: MMatrix -> IO Double
-blueNorm mm = fmap cast $ with mm c_blueNorm
-
-hypotNorm :: MMatrix -> IO Double
-hypotNorm mm = fmap cast $ with mm c_hypotNorm
-
-determinant :: MMatrix -> IO Double
-determinant mm = fmap cast $ with mm c_determinant
-
-with :: MMatrix -> (Ptr C_MatrixXd -> IO a) -> IO a
-with (MMatrix fp) f = withForeignPtr fp f
diff --git a/Data/Eigen/Matrix/Mutable.o b/Data/Eigen/Matrix/Mutable.o
new file mode 100644
Binary files /dev/null and b/Data/Eigen/Matrix/Mutable.o differ
diff --git a/Data/Eigen/Parallel.hs b/Data/Eigen/Parallel.hs
--- a/Data/Eigen/Parallel.hs
+++ b/Data/Eigen/Parallel.hs
@@ -1,5 +1,5 @@
 {- |
-Some Eigen's algorithms can exploit the multiple cores present in your hardware. To this end, it is enough to enable OpenMP on your compiler, for instance: GCC: -fopenmp ICC: -openmp MSVC: check the respective option in the build properties. You can control the number of thread that will be used using using 'setNbThreads'
+Some Eigen's algorithms can exploit the multiple cores present in your hardware. To this end, it is enough to enable OpenMP on your compiler, for instance: GCC: -fopenmp ICC: -openmp MSVC: check the respective option in the build properties. You can control the number of thread that will be used using 'setNbThreads'
 -}
 
 module Data.Eigen.Parallel where
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2013 Oleg Sidorkin
+Copyright (c) 2013-2015 Oleg Sidorkin
 
 All rights reserved.
 
diff --git a/README b/README
deleted file mode 100644
--- a/README
+++ /dev/null
@@ -1,4 +0,0 @@
-This module provides Haskell binding for Eigen C++ library.
-Some parts of documentation strings replicates original Eigen documentation.
-
-Oleg Sidorkin 
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,18 @@
+haskell-eigen
+=============
+
+This module provides Haskell binding for Eigen C++ library.
+Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.
+Eigen home page is <http://eigen.tuxfamily.org/>.
+
+Eigen is licensed under the MPL2, which is a simple weak copyleft license. Common questions about the MPL2 are answered in the official [MPL2 FAQ].
+
+Note that currently, a few features in Eigen rely on third-party code licensed under the LGPL: SimplicialCholesky, AMD ordering, and constrained\_cg. Such features are explicitly disabled by compiling Eigen with the EIGEN\_MPL2\_ONLY preprocessor symbol defined.
+
+Virtually any software may use Eigen. For example, closed-source software may use Eigen without having to disclose its own source code. Many proprietary and closed-source software projects are using Eigen right now, as well as many BSD-licensed projects.
+
+Note that **ghci** may file with *"unknown symbol ___dso_handle"* error due to dynamic linking with C++ runtime. Consider to use **ghc --make** instead. Please let me know if you know any workaround for this issue.
+
+Some parts of documentation strings are based or replicate original Eigen documentation which is available at <http://eigen.tuxfamily.org/dox/>.
+
+[MPL2 FAQ]: http://www.mozilla.org/MPL/2.0/FAQ.html
diff --git a/eigen-proxy.cpp b/eigen-proxy.cpp
--- a/eigen-proxy.cpp
+++ b/eigen-proxy.cpp
@@ -5,94 +5,153 @@
 #include <sstream>
 
 struct eigen_assert_exception : std::exception {
-	std::string _what;
-	eigen_assert_exception(const std::string& what) : _what(what) {}
-	~eigen_assert_exception() throw () {}
-	const char* what() const throw () { return _what.c_str(); }
+    std::string _what;
+    eigen_assert_exception(const std::string& what) : _what(what) {}
+    ~eigen_assert_exception() throw () {}
+    const char* what() const throw () { return _what.c_str(); }
 };
 
 void eigen_assert_fail(const char* condition, const char* function, const char* file, int line) {
-	std::ostringstream os;
+    std::ostringstream os;
     os << "assertion failed: " << condition << " in function " << function << " at " << file << ":" << line << std::endl;
     throw eigen_assert_exception(os.str());
 }
 
-MatrixXd* eigen_create(int rows, int cols) { return new MatrixXd(rows, cols); }
-MatrixXd* eigen_clone(const MatrixXd& m) { return new MatrixXd(m); }
-void eigen_destroy(const MatrixXd* m) { delete m; }
-const char* eigen_set(MatrixXd& dst, int row, int col, double val) { GUARD_START dst(row,col) = val; GUARD_END }
-const char* eigen_get(double&r, const MatrixXd& m, int row, int col) { GUARD_START r = m(row,col); GUARD_END }
-int eigen_rows(const MatrixXd& m) { return m.rows(); } 
-int eigen_cols(const MatrixXd& m) { return m.cols(); }
-double* eigen_data(MatrixXd& m) { return m.data(); } 
-void eigen_copy(MatrixXd& dst, const MatrixXd& src) { dst = src; }
-void eigen_resize(MatrixXd& m, int rows, int cols) { m.resize(rows, cols); }
-const char* eigen_add(MatrixXd& ret, const MatrixXd& lhs, const MatrixXd& rhs) { GUARD_START ret = lhs + rhs; GUARD_END }
-const char* eigen_sub(MatrixXd& ret, const MatrixXd& lhs, const MatrixXd& rhs) { GUARD_START ret = lhs - rhs; GUARD_END }
-const char* eigen_mul(MatrixXd& ret, const MatrixXd& lhs, const MatrixXd& rhs) { GUARD_START ret = lhs * rhs; GUARD_END }
-double eigen_norm(const MatrixXd& m) { return m.norm(); }
-double eigen_squaredNorm(const MatrixXd& m) { return m.squaredNorm(); }
-double eigen_blueNorm(const MatrixXd& m) { return m.squaredNorm(); }
-// double mstableNorm(const MatrixXd& m) { return m.stableNorm(); }
-double eigen_hypotNorm(const MatrixXd& m) { return m.hypotNorm(); }
-double eigen_sum(const MatrixXd& m) { return m.sum(); }
-double eigen_prod(const MatrixXd& m) { return m.prod(); }
-double eigen_mean(const MatrixXd& m) { return m.mean(); }
-double eigen_minCoeff(const MatrixXd& m) { return m.minCoeff(); }
-double eigen_maxCoeff(const MatrixXd& m) { return m.maxCoeff(); }
-double eigen_trace(const MatrixXd& m) { return m.trace(); }
+typedef Map< Matrix<double,Dynamic,Dynamic> > MapMatrix;
 
-const char* eigen_inverse(MatrixXd& m) { GUARD_START m = m.inverse().eval(); GUARD_END }
-const char* eigen_adjoint(MatrixXd& m) { GUARD_START m.adjointInPlace(); GUARD_END }
-const char* eigen_transpose(MatrixXd& m) { GUARD_START m.transposeInPlace(); GUARD_END }
-const char* eigen_normalize(MatrixXd& m) { GUARD_START m.normalize(); GUARD_END }
-const char* eigen_conjugate(MatrixXd& m) { GUARD_START m = m.conjugate().eval(); GUARD_END }
+extern "C" {
 
+const char* eigen_add(double* data, int rows, int cols, const double* data1, int rows1, int cols1) 
+{
+    GUARD_START
+    MapMatrix(data, rows, cols) += MapMatrix(data1, rows1, cols1);
+    GUARD_END
+}
 
-double eigen_determinant(const MatrixXd& m) { return m.determinant(); }
+const char* eigen_sub(double* data, int rows, int cols, const double* data1, int rows1, int cols1) 
+{
+    GUARD_START
+    MapMatrix(data, rows, cols) -= MapMatrix(data1, rows1, cols1);
+    GUARD_END
+}
 
-const char* eigen_solve(Decomposition d, MatrixXd& x, const MatrixXd& A, const MatrixXd& b) {
-	GUARD_START
-	switch (d) {
-		case ::PartialPivLU:
-			x = A.partialPivLu().solve(b); 
-			break;
-		case ::FullPivLU:
-			x = A.fullPivLu().solve(b); 
-			break;
-		case ::HouseholderQR:
-			x = A.householderQr().solve(b); 
-			break;
-		case ::ColPivHouseholderQR:
-			x = A.colPivHouseholderQr().solve(b); 
-			break;
-		case ::FullPivHouseholderQR: 
-			x = A.fullPivHouseholderQr().solve(b); 
-			break;
-		case ::LLT:
-			x = A.llt().solve(b); 
-			break;
-		case ::LDLT:
-			x = A.ldlt().solve(b); 
-			break;
-		case ::JacobiSVD:
-			x = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b); 
-			break;
-	}
-	GUARD_END
+const char* eigen_mul(double* data, int rows, int cols, const double* data1, int rows1, int cols1)  
+{
+    GUARD_START
+    MapMatrix(data, rows, cols) *= MapMatrix(data1, rows1, cols1);
+    GUARD_END
 }
 
-const char* eigen_relativeError(double& e, const MatrixXd& x, const MatrixXd& A, const MatrixXd& b) {
-	GUARD_START
-	e = (A*x - b).norm() / b.norm();
-	GUARD_END
+double eigen_norm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).norm();  }
+double eigen_squaredNorm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).squaredNorm(); }
+double eigen_blueNorm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).squaredNorm(); }
+double eigen_hypotNorm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).hypotNorm(); }
+double eigen_sum(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).sum(); }
+double eigen_prod(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).prod(); }
+double eigen_mean(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).mean(); }
+double eigen_trace(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).trace(); }
+double eigen_determinant(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).determinant(); }
+
+const char* eigen_inverse(double* data, int rows, int cols, const double* data1, int rows1, int cols1) 
+{
+    GUARD_START
+    puts("inverse");
+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).inverse();
+    GUARD_END
 }
 
+const char* eigen_adjoint(double* data, int rows, int cols, const double* data1, int rows1, int cols1) 
+{
+    GUARD_START
+    puts("adjoint");
+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).adjoint();
+    GUARD_END
+}
+
+const char* eigen_conjugate(double* data, int rows, int cols, const double* data1, int rows1, int cols1) 
+{
+    GUARD_START
+    puts("conjugate");
+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).conjugate();
+    GUARD_END
+}
+
+const char* eigen_transpose(double* data, int rows, int cols, const double* data1, int rows1, int cols1) 
+{
+    GUARD_START
+    puts("transpose");
+    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).transpose();
+    GUARD_END
+}
+
+const char* eigen_normalize(double* data, int rows, int cols) 
+{
+    GUARD_START
+    puts("normalize");
+    MapMatrix(data, rows, cols).normalize();
+    GUARD_END
+}
+
+
+const char* eigen_solve(Decomposition d, 
+    double* px, int rx, int cx, // x
+    const double* pa, int ra, int ca, // A
+    const double* pb, int rb, int cb) // b
+{
+    GUARD_START
+    MapMatrix x(px, rx, cx);
+    MapMatrix A(pa, ra, ca);
+    MapMatrix b(pb, rb, cb);
+    switch (d) {
+        case ::PartialPivLU:
+            x = A.partialPivLu().solve(b);
+            break;
+        case ::FullPivLU:
+            x = A.fullPivLu().solve(b);
+            break;
+        case ::HouseholderQR:
+            x = A.householderQr().solve(b);
+            break;
+        case ::ColPivHouseholderQR:
+            x = A.colPivHouseholderQr().solve(b);
+            break;
+        case ::FullPivHouseholderQR:
+            x = A.fullPivHouseholderQr().solve(b);
+            break;
+        case ::LLT:
+            x = A.llt().solve(b);
+            break;
+        case ::LDLT:
+            x = A.ldlt().solve(b);
+            break;
+        case ::JacobiSVD:
+            x = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b);
+            break;
+    }
+    GUARD_END
+}
+
+
+const char* eigen_relativeError(double& e, 
+    const double* px, int rx, int cx, // x
+    const double* pa, int ra, int ca, // A
+    const double* pb, int rb, int cb) // b
+{
+    GUARD_START
+    MapMatrix x(px, rx, cx);
+    MapMatrix A(pa, ra, ca);
+    MapMatrix b(pb, rb, cb);
+    e = (A*x - b).norm() / b.norm();
+    GUARD_END
+}
+
+
 void eigen_initParallel() {
-	initParallel();
+    initParallel();
 }
 
 void eigen_setNbThreads(int n) {
-	setNbThreads(n);
+    setNbThreads(n);
 }
 
+}
diff --git a/eigen-proxy.h b/eigen-proxy.h
--- a/eigen-proxy.h
+++ b/eigen-proxy.h
@@ -11,35 +11,23 @@
 
 extern "C" {
 
-MatrixXd* eigen_create(int rows, int cols);
-void eigen_destroy(const MatrixXd*);
-MatrixXd* eigen_clone(const MatrixXd& m);
-const char* eigen_set(MatrixXd&,int,int,double);
-const char* eigen_get(double&,const MatrixXd&,int,int);
-int eigen_rows(const MatrixXd&);
-int eigen_cols(const MatrixXd&);
-double* eigen_data(MatrixXd&);
-void eigen_copy(MatrixXd&, const MatrixXd&);
-void eigen_resize(MatrixXd&, int, int);
-const char* eigen_add(MatrixXd&, const MatrixXd&, const MatrixXd&);
-const char* eigen_sub(MatrixXd&, const MatrixXd&, const MatrixXd&);
-const char* eigen_mul(MatrixXd&, const MatrixXd&, const MatrixXd&);
-double eigen_norm(const MatrixXd&);
-double eigen_squaredNorm(const MatrixXd&);
-double eigen_blueNorm(const MatrixXd&);
-double eigen_stableNorm(const MatrixXd&);
-double eigen_hypotNorm(const MatrixXd&);
-double eigen_sum(const MatrixXd&);
-double eigen_prod(const MatrixXd&);
-double eigen_mean(const MatrixXd&);
-double eigen_minCoeff(const MatrixXd&);
-double eigen_maxCoeff(const MatrixXd&);
-double eigen_trace(const MatrixXd&);
-const char* eigen_transpose(MatrixXd&);
-const char* eigen_normalize(MatrixXd&);
-const char* eigen_inverse(MatrixXd&);
-const char* eigen_adjoint(MatrixXd&);
-double eigen_determinant(const MatrixXd& src);
+const char* eigen_add(double*, int, int, const double*, int, int);
+const char* eigen_sub(double*, int, int, const double*, int, int);
+const char* eigen_mul(double*, int, int, const double*, int, int);
+double eigen_norm(const double*, int, int);
+double eigen_squaredNorm(const double*, int, int);
+double eigen_blueNorm(const double*, int, int);
+double eigen_hypotNorm(const double*, int, int);
+double eigen_sum(const double*, int, int);
+double eigen_prod(const double*, int, int);
+double eigen_mean(const double*, int, int);
+double eigen_trace(const double*, int, int);
+double eigen_determinant( const double*, int, int);
+const char* eigen_transpose(double*, int, int, const double*, int, int);
+const char* eigen_normalize(double*, int, int);
+const char* eigen_inverse(double*, int, int, const double*, int, int);
+const char* eigen_conjugate(double*, int, int, const double*, int, int);
+const char* eigen_adjoint(double*, int, int, const double*, int, int);
 void eigen_initParallel();
 void eigen_setNbThreads(int);
 
@@ -47,8 +35,15 @@
 	PartialPivLU, FullPivLU, HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR, LLT, LDLT, JacobiSVD
 };
 
-const char* eigen_solve(Decomposition d, MatrixXd& x, const MatrixXd& A, const MatrixXd& b); // Ax=b
-const char* eigen_relativeError(double& e, const MatrixXd& x, const MatrixXd& A, const MatrixXd& b);
+const char* eigen_solve(Decomposition d, // Ax=b 
+	double*, int, int, // x
+	const double*, int, int, // A
+	const double*, int, int); // b
+
+const char* eigen_relativeError(double& e, 
+	const double*, int, int, // x
+	const double*, int, int, // A
+	const double*, int, int); // b;
 
 
 } // end extern "C"
diff --git a/eigen-proxy.o b/eigen-proxy.o
new file mode 100644
# file too large to diff: eigen-proxy.o
diff --git a/eigen.cabal b/eigen.cabal
--- a/eigen.cabal
+++ b/eigen.cabal
@@ -1,7 +1,7 @@
 name:           eigen
-version:        1.0.0
+version:        1.1.1
 homepage:       https://github.com/osidorkin/haskell-eigen
-synopsis:       Haskel binding for Eigen library
+synopsis:       Eigen C++ library (linear algebra: matrices, vectors, numerical solvers, and related algorithms).
 description:    This module provides Haskell binding for Eigen C++ library.
                 Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.
                 Eigen home page is <http://eigen.tuxfamily.org/>.
@@ -12,10 +12,8 @@
                 .
                 Virtually any software may use Eigen. For example, closed-source software may use Eigen without having to disclose its own source code. Many proprietary and closed-source software projects are using Eigen right now, as well as many BSD-licensed projects.
                 .
-                Note that ghci may fail with \"@unknown symbol ___dso_handle@\" error due to dynamic linking with C++ runtime. Consider to use ghc --make instead. Please let me know if you know any workaround for this issue.
-                .
                 Some parts of documentation strings are based or replicate original Eigen documentation which is available at <http://eigen.tuxfamily.org/dox/>.
-category:       Data
+category:       Data, Math, Algebra, Statistics, Algorithms
 license:        BSD3
 license-file:   LICENSE
 copyright:      (c) 2013, Oleg Sidorkin
@@ -23,7 +21,7 @@
 maintainer:     Oleg Sidorkin <oleg.sidorkin@gmail.com>
 cabal-version:  >= 1.6
 build-type:     Simple
-tested-with:    GHC ==7.4.1
+tested-with:    GHC ==7.8.3
 
 source-repository head
   type:     git
@@ -39,7 +37,8 @@
 
     ghc-options:        -Wall -fno-warn-name-shadowing
     build-depends:      base >= 3 && < 5,
-                        vector >= 0.5 && < 0.11
+                        vector >= 0.5 && < 0.11,
+                        primitive
 
     include-dirs:       eigen3
     c-sources:          eigen-proxy.cpp
diff --git a/eigen3/Eigen/README b/eigen3/Eigen/README
deleted file mode 100644
--- a/eigen3/Eigen/README
+++ /dev/null
@@ -1,4 +0,0 @@
-This folder contains headers and actual licenses from eigen-3.1.2 distribution.
-The original distribution package is available at http://bitbucket.org/eigen/eigen/get/3.1.2.tar.bz2
-
-Oleg Sidorkin
diff --git a/eigen3/README b/eigen3/README
new file mode 100644
--- /dev/null
+++ b/eigen3/README
@@ -0,0 +1,4 @@
+This folder contains headers and actual licenses from eigen-3.1.2 distribution.
+The original distribution package is available at http://bitbucket.org/eigen/eigen/get/3.1.2.tar.bz2
+
+Oleg Sidorkin