diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) 2009, 2010 Bryan O'Sullivan
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,5 @@
+# statistics-dense-linear-algebra
+
+[![Build Status](https://travis-ci.org/githubuser/statistics-dense-linear-algebra.png)](https://travis-ci.org/githubuser/statistics-dense-linear-algebra)
+
+The dense linear algebra functionality and related modules, extracted from `statistics-0.14.0.2`
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/dense-linear-algebra.cabal b/dense-linear-algebra.cabal
new file mode 100644
--- /dev/null
+++ b/dense-linear-algebra.cabal
@@ -0,0 +1,50 @@
+name:                dense-linear-algebra
+version:             0.1.0.0
+synopsis:            Simple and incomplete pure haskell implementation of linear algebra
+description:
+  This library is simply collection of linear-algebra related modules
+  split from statistics library.
+
+license:             BSD2
+license-file:        LICENSE
+author:              Bryan O'Sullivan <bos@serpentine.com>
+maintainer:          Alexey Khudaykov <alexey.skladnoy@gmail.com>
+copyright:           2018 Author name here
+category:            Math, Statistics, Numeric
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+tested-with: GHC==7.6.3, GHC==7.8.3, GHC==7.10.3, GHC==8.0.1
+
+library
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  hs-source-dirs:      src
+  exposed-modules:     Statistics.Matrix
+                       Statistics.Matrix.Algorithms
+                       Statistics.Matrix.Function
+                       Statistics.Matrix.Mutable
+                       Statistics.Matrix.Types
+  build-depends:       base                    >= 4.5 && < 5
+                     , deepseq                 >= 1.1.0.2
+                     , math-functions          >= 0.1.7                     
+                     , primitive               >= 0.3
+                     , vector                  >= 0.10
+                     , vector-algorithms       >= 0.4
+                     , vector-th-unbox
+                     , vector-binary-instances >= 0.2.1
+
+test-suite spec
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             LibSpec.hs
+  build-depends:       base
+                     , dense-linear-algebra
+                     , hspec
+                     , QuickCheck
+
+source-repository head
+  type:     git
+  location: https://github.com/bos/statistics
diff --git a/src/Statistics/Matrix.hs b/src/Statistics/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Matrix.hs
@@ -0,0 +1,265 @@
+{-# LANGUAGE PatternGuards #-}
+-- |
+-- Module    : Statistics.Matrix
+-- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan
+-- License   : BSD3
+--
+-- Basic matrix operations.
+--
+-- There isn't a widely used matrix package for Haskell yet, so
+-- we implement the necessary minimum here.
+
+module Statistics.Matrix
+    ( -- * Data types
+      Matrix(..)
+    , Vector
+      -- * Conversion from/to lists/vectors
+    , fromVector
+    , fromList
+    , fromRowLists
+    , fromRows
+    , fromColumns
+    , toVector
+    , toList
+    , toRows
+    , toColumns
+    , toRowLists
+      -- * Other
+    , generate
+    , generateSym
+    , ident
+    , diag
+    , dimension
+    , center
+    , multiply
+    , multiplyV
+    , transpose
+    , power
+    , norm
+    , column
+    , row
+    , map
+    , for
+    , unsafeIndex
+    , hasNaN
+    , bounds
+    , unsafeBounds
+    ) where
+
+import Prelude hiding (exponent, map)
+import Control.Applicative ((<$>))
+import Control.Monad.ST
+import qualified Data.Vector.Unboxed as U
+import           Data.Vector.Unboxed   ((!))
+import qualified Data.Vector.Unboxed.Mutable as UM
+
+import Numeric.Sum         (sumVector,kbn)
+import Statistics.Matrix.Function
+import Statistics.Matrix.Types
+import Statistics.Matrix.Mutable  (unsafeNew,unsafeWrite,unsafeFreeze)
+
+
+
+----------------------------------------------------------------
+-- Conversion to/from vectors/lists
+----------------------------------------------------------------
+
+-- | Convert from a row-major list.
+fromList :: Int                 -- ^ Number of rows.
+         -> Int                 -- ^ Number of columns.
+         -> [Double]            -- ^ Flat list of values, in row-major order.
+         -> Matrix
+fromList r c = fromVector r c . U.fromList
+
+-- | create a matrix from a list of lists, as rows
+fromRowLists :: [[Double]] -> Matrix
+fromRowLists = fromRows . fmap U.fromList
+
+-- | Convert from a row-major vector.
+fromVector :: Int               -- ^ Number of rows.
+           -> Int               -- ^ Number of columns.
+           -> U.Vector Double   -- ^ Flat list of values, in row-major order.
+           -> Matrix
+fromVector r c v
+  | r*c /= len = error "input size mismatch"
+  | otherwise  = Matrix r c v
+  where len    = U.length v
+
+-- | create a matrix from a list of vectors, as rows
+fromRows :: [Vector] -> Matrix
+fromRows xs
+  | [] <- xs        = error "Statistics.Matrix.fromRows: empty list of rows!"
+  | any (/=nCol) ns = error "Statistics.Matrix.fromRows: row sizes do not match"
+  | nCol == 0       = error "Statistics.Matrix.fromRows: zero columns in matrix"
+  | otherwise       = fromVector nRow nCol (U.concat xs)
+  where
+    nCol:ns = U.length <$> xs
+    nRow    = length xs
+
+
+-- | create a matrix from a list of vectors, as columns
+fromColumns :: [Vector] -> Matrix
+fromColumns = transpose . fromRows
+
+-- | Convert to a row-major flat vector.
+toVector :: Matrix -> U.Vector Double
+toVector (Matrix _ _ v) = v
+
+-- | Convert to a row-major flat list.
+toList :: Matrix -> [Double]
+toList = U.toList . toVector
+
+-- | Convert to a list of lists, as rows
+toRowLists :: Matrix -> [[Double]]
+toRowLists (Matrix _ nCol v)
+  = chunks $ U.toList v
+  where
+    chunks [] = []
+    chunks xs = case splitAt nCol xs of
+      (rowE,rest) -> rowE : chunks rest
+
+
+-- | Convert to a list of vectors, as rows
+toRows :: Matrix -> [Vector]
+toRows (Matrix _ nCol v) = chunks v
+  where
+    chunks xs
+      | U.null xs = []
+      | otherwise = case U.splitAt nCol xs of
+          (rowE,rest) -> rowE : chunks rest
+
+-- | Convert to a list of vectors, as columns
+toColumns :: Matrix -> [Vector]
+toColumns = toRows . transpose
+
+
+
+----------------------------------------------------------------
+-- Other
+----------------------------------------------------------------
+
+-- | Generate matrix using function
+generate :: Int                 -- ^ Number of rows
+         -> Int                 -- ^ Number of columns
+         -> (Int -> Int -> Double)
+            -- ^ Function which takes /row/ and /column/ as argument.
+         -> Matrix
+generate nRow nCol f
+  = Matrix nRow nCol $ U.generate (nRow*nCol) $ \i ->
+      let (r,c) = i `quotRem` nCol in f r c
+
+-- | Generate symmetric square matrix using function
+generateSym
+  :: Int                 -- ^ Number of rows and columns
+  -> (Int -> Int -> Double)
+     -- ^ Function which takes /row/ and /column/ as argument. It must
+     --   be symmetric in arguments: @f i j == f j i@
+  -> Matrix
+generateSym n f = runST $ do
+  m <- unsafeNew n n
+  for 0 n $ \r -> do
+    unsafeWrite m r r (f r r)
+    for (r+1) n $ \c -> do
+      let x = f r c
+      unsafeWrite m r c x
+      unsafeWrite m c r x
+  unsafeFreeze m
+
+
+-- | Create the square identity matrix with given dimensions.
+ident :: Int -> Matrix
+ident n = diag $ U.replicate n 1.0
+
+-- | Create a square matrix with given diagonal, other entries default to 0
+diag :: Vector -> Matrix
+diag v
+  = Matrix n n $ U.create $ do
+      arr <- UM.replicate (n*n) 0
+      for 0 n $ \i ->
+        UM.unsafeWrite arr (i*n + i) (v ! i)
+      return arr
+  where
+    n = U.length v
+
+-- | Return the dimensions of this matrix, as a (row,column) pair.
+dimension :: Matrix -> (Int, Int)
+dimension (Matrix r c _) = (r, c)
+
+-- | Matrix-matrix multiplication. Matrices must be of compatible
+-- sizes (/note: not checked/).
+multiply :: Matrix -> Matrix -> Matrix
+multiply m1@(Matrix r1 _ _) m2@(Matrix _ c2 _) =
+  Matrix r1 c2 $ U.generate (r1*c2) go
+  where
+    go t = sumVector kbn $ U.zipWith (*) (row m1 i) (column m2 j)
+      where (i,j) = t `quotRem` c2
+
+-- | Matrix-vector multiplication.
+multiplyV :: Matrix -> Vector -> Vector
+multiplyV m v
+  | cols m == c = U.generate (rows m) (sumVector kbn . U.zipWith (*) v . row m)
+  | otherwise   = error $ "matrix/vector unconformable " ++ show (cols m,c)
+  where c = U.length v
+
+-- | Raise matrix to /n/th power. Power must be positive
+-- (/note: not checked).
+power :: Matrix -> Int -> Matrix
+power mat 1 = mat
+power mat n = res
+  where
+    mat2 = power mat (n `quot` 2)
+    pow  = multiply mat2 mat2
+    res | odd n     = multiply pow mat
+        | otherwise = pow
+
+-- | Element in the center of matrix (not corrected for exponent).
+center :: Matrix -> Double
+center mat@(Matrix r c _) =
+    unsafeBounds U.unsafeIndex mat (r `quot` 2) (c `quot` 2)
+
+-- | Calculate the Euclidean norm of a vector.
+norm :: Vector -> Double
+norm = sqrt . sumVector kbn . U.map square
+
+-- | Return the given column.
+column :: Matrix -> Int -> Vector
+column (Matrix r c v) i = U.backpermute v $ U.enumFromStepN i c r
+{-# INLINE column #-}
+
+-- | Return the given row.
+row :: Matrix -> Int -> Vector
+row (Matrix _ c v) i = U.slice (c*i) c v
+
+unsafeIndex :: Matrix
+            -> Int              -- ^ Row.
+            -> Int              -- ^ Column.
+            -> Double
+unsafeIndex = unsafeBounds U.unsafeIndex
+
+-- | Apply function to every element of matrix
+map :: (Double -> Double) -> Matrix -> Matrix
+map f (Matrix r c v) = Matrix r c (U.map f v)
+
+-- | Indicate whether any element of the matrix is @NaN@.
+hasNaN :: Matrix -> Bool
+hasNaN = U.any isNaN . toVector
+
+-- | Given row and column numbers, calculate the offset into the flat
+-- row-major vector.
+bounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
+bounds k (Matrix rs cs v) r c
+  | r < 0 || r >= rs = error "row out of bounds"
+  | c < 0 || c >= cs = error "column out of bounds"
+  | otherwise        = k v $! r * cs + c
+{-# INLINE bounds #-}
+
+-- | Given row and column numbers, calculate the offset into the flat
+-- row-major vector, without checking.
+unsafeBounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
+unsafeBounds k (Matrix _ cs v) r c = k v $! r * cs + c
+{-# INLINE unsafeBounds #-}
+
+transpose :: Matrix -> Matrix
+transpose m@(Matrix r0 c0 _) = Matrix c0 r0 . U.generate (r0*c0) $ \i ->
+  let (r,c) = i `quotRem` r0
+  in unsafeIndex m c r
diff --git a/src/Statistics/Matrix/Algorithms.hs b/src/Statistics/Matrix/Algorithms.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Matrix/Algorithms.hs
@@ -0,0 +1,42 @@
+-- |
+-- Module    : Statistics.Matrix.Algorithms
+-- Copyright : 2014 Bryan O'Sullivan
+-- License   : BSD3
+--
+-- Useful matrix functions.
+
+module Statistics.Matrix.Algorithms
+    (
+      qr
+    ) where
+
+import Control.Applicative ((<$>), (<*>))
+import Control.Monad.ST (ST, runST)
+import Prelude hiding (replicate)
+import Numeric.Sum       (sumVector,kbn)
+import Statistics.Matrix (Matrix, column, dimension, for, norm)
+import qualified Statistics.Matrix.Mutable as M
+import qualified Data.Vector.Unboxed as U
+
+-- | /O(r*c)/ Compute the QR decomposition of a matrix.
+-- The result returned is the matrices (/q/,/r/).
+qr :: Matrix -> (Matrix, Matrix)
+qr mat = runST $ do
+  let (m,n) = dimension mat
+  r <- M.replicate n n 0
+  a <- M.thaw mat
+  for 0 n $ \j -> do
+    cn <- M.immutably a $ \aa -> norm (column aa j)
+    M.unsafeWrite r j j cn
+    for 0 m $ \i -> M.unsafeModify a i j (/ cn)
+    for (j+1) n $ \jj -> do
+      p <- innerProduct a j jj
+      M.unsafeWrite r j jj p
+      for 0 m $ \i -> do
+        aij <- M.unsafeRead a i j
+        M.unsafeModify a i jj $ subtract (p * aij)
+  (,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r
+
+innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double
+innerProduct mmat j k = M.immutably mmat $ \mat ->
+  sumVector kbn $ U.zipWith (*) (column mat j) (column mat k)
diff --git a/src/Statistics/Matrix/Function.hs b/src/Statistics/Matrix/Function.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Matrix/Function.hs
@@ -0,0 +1,15 @@
+{-# LANGUAGE BangPatterns #-}
+-- |
+module Statistics.Matrix.Function where
+
+-- | Multiply a number by itself.
+square :: Double -> Double
+square x = x * x
+
+-- | Simple for loop.  Counts from /start/ to /end/-1.
+for :: Monad m => Int -> Int -> (Int -> m ()) -> m ()
+for n0 !n f = loop n0
+  where
+    loop i | i == n    = return ()
+           | otherwise = f i >> loop (i+1)
+{-# INLINE for #-}
diff --git a/src/Statistics/Matrix/Mutable.hs b/src/Statistics/Matrix/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Matrix/Mutable.hs
@@ -0,0 +1,86 @@
+-- |
+-- Module    : Statistics.Matrix.Mutable
+-- Copyright : (c) 2014 Bryan O'Sullivan
+-- License   : BSD3
+--
+-- Basic mutable matrix operations.
+
+module Statistics.Matrix.Mutable
+    (
+      MMatrix(..)
+    , MVector
+    , replicate
+    , thaw
+    , bounds
+    , unsafeNew
+    , unsafeFreeze
+    , unsafeRead
+    , unsafeWrite
+    , unsafeModify
+    , immutably
+    , unsafeBounds
+    ) where
+
+import Control.Applicative ((<$>))
+import Control.DeepSeq (NFData(..))
+import Control.Monad.ST (ST)
+import Statistics.Matrix.Types (Matrix(..), MMatrix(..), MVector)
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as M
+import Prelude hiding (replicate)
+
+replicate :: Int -> Int -> Double -> ST s (MMatrix s)
+replicate r c k = MMatrix r c <$> M.replicate (r*c) k
+
+thaw :: Matrix -> ST s (MMatrix s)
+thaw (Matrix r c v) = MMatrix r c <$> U.thaw v
+
+unsafeFreeze :: MMatrix s -> ST s Matrix
+unsafeFreeze (MMatrix r c mv) = Matrix r c <$> U.unsafeFreeze mv
+
+-- | Allocate new matrix. Matrix content is not initialized hence unsafe.
+unsafeNew :: Int                -- ^ Number of row
+          -> Int                -- ^ Number of columns
+          -> ST s (MMatrix s)
+unsafeNew r c
+  | r < 0     = error "Statistics.Matrix.Mutable.unsafeNew: negative number of rows"
+  | c < 0     = error "Statistics.Matrix.Mutable.unsafeNew: negative number of columns"
+  | otherwise = do
+      vec <- M.new (r*c)
+      return $ MMatrix r c vec
+
+unsafeRead :: MMatrix s -> Int -> Int -> ST s Double
+unsafeRead mat r c = unsafeBounds mat r c M.unsafeRead
+{-# INLINE unsafeRead #-}
+
+unsafeWrite :: MMatrix s -> Int -> Int -> Double -> ST s ()
+unsafeWrite mat row col k = unsafeBounds mat row col $ \v i ->
+  M.unsafeWrite v i k
+{-# INLINE unsafeWrite #-}
+
+unsafeModify :: MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
+unsafeModify mat row col f = unsafeBounds mat row col $ \v i -> do
+  k <- M.unsafeRead v i
+  M.unsafeWrite v i (f k)
+{-# INLINE unsafeModify #-}
+
+-- | Given row and column numbers, calculate the offset into the flat
+-- row-major vector.
+bounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r
+bounds (MMatrix rs cs mv) r c k
+  | r < 0 || r >= rs = error "row out of bounds"
+  | c < 0 || c >= cs = error "column out of bounds"
+  | otherwise        = k mv $! r * cs + c
+{-# INLINE bounds #-}
+
+-- | Given row and column numbers, calculate the offset into the flat
+-- row-major vector, without checking.
+unsafeBounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r
+unsafeBounds (MMatrix _ cs mv) r c k = k mv $! r * cs + c
+{-# INLINE unsafeBounds #-}
+
+immutably :: NFData a => MMatrix s -> (Matrix -> a) -> ST s a
+immutably mmat f = do
+  k <- f <$> unsafeFreeze mmat
+  rnf k `seq` return k
+{-# INLINE immutably #-}
diff --git a/src/Statistics/Matrix/Types.hs b/src/Statistics/Matrix/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Statistics/Matrix/Types.hs
@@ -0,0 +1,60 @@
+-- |
+-- Module    : Statistics.Matrix.Types
+-- Copyright : 2014 Bryan O'Sullivan
+-- License   : BSD3
+--
+-- Basic matrix operations.
+--
+-- There isn't a widely used matrix package for Haskell yet, so
+-- we implement the necessary minimum here.
+
+module Statistics.Matrix.Types
+    (
+      Vector
+    , MVector
+    , Matrix(..)
+    , MMatrix(..)
+    , debug
+    ) where
+
+import Data.Char (isSpace)
+import Numeric (showFFloat)
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as M
+
+type Vector = U.Vector Double
+type MVector s = M.MVector s Double
+
+-- | Two-dimensional matrix, stored in row-major order.
+data Matrix = Matrix {
+      rows     :: {-# UNPACK #-} !Int -- ^ Rows of matrix.
+    , cols     :: {-# UNPACK #-} !Int -- ^ Columns of matrix.
+    , _vector  :: !Vector             -- ^ Matrix data.
+    } deriving (Eq)
+
+-- | Two-dimensional mutable matrix, stored in row-major order.
+data MMatrix s = MMatrix
+                 {-# UNPACK #-} !Int
+                 {-# UNPACK #-} !Int
+                 !(MVector s)
+
+-- The Show instance is useful only for debugging.
+instance Show Matrix where
+    show = debug
+
+debug :: Matrix -> String
+debug (Matrix r c vs) = unlines $ zipWith (++) (hdr0 : repeat hdr) rrows
+  where
+    rrows         = map (cleanEnd . unwords) . split $ zipWith (++) ldone tdone
+    hdr0          = show (r,c) ++ " "
+    hdr           = replicate (length hdr0) ' '
+    pad plus k xs = replicate (k - length xs) ' ' `plus` xs
+    ldone         = map (pad (++) (longest lstr)) lstr
+    tdone         = map (pad (flip (++)) (longest tstr)) tstr
+    (lstr, tstr)  = unzip . map (break (=='.') . render) . U.toList $ vs
+    longest       = maximum . map length
+    render k      = reverse . dropWhile (=='.') . dropWhile (=='0') . reverse .
+                    showFFloat (Just 4) k $ ""
+    split []      = []
+    split xs      = i : split rest where (i, rest) = splitAt c xs
+    cleanEnd      = reverse . dropWhile isSpace . reverse
diff --git a/test/LibSpec.hs b/test/LibSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/LibSpec.hs
@@ -0,0 +1,16 @@
+module Main where
+
+import Test.Hspec
+import Test.Hspec.QuickCheck
+
+
+main :: IO ()
+main = hspec spec
+
+spec :: Spec
+spec =
+  describe "Lib" $ do
+    it "works" $ do
+      True `shouldBe` True
+    -- prop "ourAdd is commutative" $ \x y ->
+    --   ourAdd x y `shouldBe` ourAdd y x