diff --git a/Data/NeuralNetwork/Backend/HMatrix.hs b/Data/NeuralNetwork/Backend/HMatrix.hs
new file mode 100644
--- /dev/null
+++ b/Data/NeuralNetwork/Backend/HMatrix.hs
@@ -0,0 +1,91 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeFamilies #-}
+module Data.NeuralNetwork.Backend.HMatrix (
+  module Data.NeuralNetwork.Backend.HMatrix.Layers,
+  ByHmatrix(..),
+  ErrCode(..)
+) where
+
+import Data.NeuralNetwork
+import Data.NeuralNetwork.Backend.HMatrix.Utils
+import Data.NeuralNetwork.Backend.HMatrix.Layers
+import Numeric.LinearAlgebra (Vector, Matrix)
+import Control.Monad.Except
+import Data.Functor.Identity
+
+data ErrCode = ErrMismatch
+type Err     = ExceptT ErrCode IO
+
+-- the backend type
+data ByHmatrix = ByHmatrix
+
+-- with 1D input
+instance (TranslateBody s, Component (RunLayer (SpecToTag s))) =>
+    Backend ByHmatrix (SpecIn1D :++ s) where
+  type Env ByHmatrix = Err
+  type ConvertFromSpec (SpecIn1D :++ s) = RunLayer (SpecToTag s)
+  compile _ (a :++ l)= trans (size Nothing a) l
+
+-- with 2D input
+instance (TranslateBody s, Component (RunLayer (SpecToTag s))) =>
+    Backend ByHmatrix (SpecIn2D :++ s) where
+  type Env ByHmatrix = Err
+  type ConvertFromSpec (SpecIn2D :++ s) = RunLayer (SpecToTag s)
+  compile _ (a :++ l)= trans (size Nothing a) l
+
+instance RunInEnv Identity Err where
+  run = return . runIdentity
+
+-- It is necessary to propagate the size along the layers,
+-- because fullconnect and convolution need to know
+-- the previous size.
+data Size = D1 Int | D2 Int Int Int
+
+class ComputeSize l where
+  size :: Maybe Size -> l -> Size
+instance ComputeSize SpecIn1D where
+  size Nothing (In1D n) = D1 n
+instance ComputeSize SpecIn2D where
+  size Nothing (In2D m n) = D2 1 m n
+instance ComputeSize SpecReshape2DAs1D where
+  size (Just (D2 k m n)) _ = D1 (k*m*n)
+instance ComputeSize SpecFullConnect where
+  size _ (FullConnect n)   = D1 n
+instance ComputeSize SpecConvolution where
+  size (Just (D2 _ m n)) (Convolution k f p) = D2 k (m+2*p-f+1) (n+2*p-f+1)
+instance ComputeSize SpecMaxPooling where
+  size (Just (D2 k m n)) (MaxPooling s) = D2 k (m `div` s) (n `div` s)
+
+-- translate the body of specification
+class TranslateBody s where
+  type SpecToTag s
+  trans :: Size -> s -> Err (RunLayer (SpecToTag s))
+
+instance TranslateBody SpecFullConnect where
+  type SpecToTag SpecFullConnect = S F (T (SinglC :. Vector))
+  trans (D1 s) (FullConnect n) = do u <- lift $ newFLayer s n
+                                    return $ Stack u (Activation (relu, relu'))
+  trans _ _ = throwError ErrMismatch
+
+instance TranslateBody SpecConvolution where
+  type SpecToTag SpecConvolution = S C (T (MultiC :. Matrix))
+  trans (D2 k s t) (Convolution n f p) = do u <- lift $ newCLayer k n f p
+                                            return $ Stack u (Activation (relu, relu'))
+  trans _ _ = throwError ErrMismatch
+
+instance TranslateBody SpecReshape2DAs1D where
+  type SpecToTag SpecReshape2DAs1D = A
+  trans (D2 _ _ _) _ = return As1D
+  trans (D1 _)     _ = throwError ErrMismatch
+
+instance TranslateBody SpecMaxPooling where
+  type SpecToTag SpecMaxPooling = M
+  trans (D2 _ _ _) (MaxPooling n) = return (MaxP n)
+  trans (D1 _)     _              = throwError ErrMismatch
+
+instance (TranslateBody a, TranslateBody c, ComputeSize a) => TranslateBody (a :++ c) where
+  type SpecToTag (a :++ b) = S (SpecToTag a) (SpecToTag b)
+  trans s (a :++ c) = do u <- trans s a
+                         v <- trans (size (Just s) a) c
+                         return $ Stack u v
diff --git a/Data/NeuralNetwork/Backend/HMatrix/Layers.hs b/Data/NeuralNetwork/Backend/HMatrix/Layers.hs
new file mode 100644
--- /dev/null
+++ b/Data/NeuralNetwork/Backend/HMatrix/Layers.hs
@@ -0,0 +1,345 @@
+{-# LANGUAGE BangPatterns, TypeFamilies, TypeOperators, FlexibleInstances, FlexibleContexts, GADTs #-}
+module Data.NeuralNetwork.Backend.HMatrix.Layers where
+
+import Numeric.LinearAlgebra hiding (R, C)
+import Numeric.LinearAlgebra.Devel
+import qualified Data.Vector as V
+import qualified Data.Vector.Storable as SV
+import qualified Data.Vector.Storable.Mutable as SVM
+import System.Random.MWC
+import System.Random.MWC.Distributions
+import Control.Monad.ST
+import Control.Monad (liftM2, forM_, when)
+import GHC.Float
+import Data.STRef
+import Data.Functor.Identity
+import Control.DeepSeq
+import Data.NeuralNetwork
+import Data.NeuralNetwork.Backend.HMatrix.Utils
+
+type R = Float
+
+-- We parameterise the activation layer T, where the parameter indicates how
+-- elements are contained:
+--   SinglC :. Vector, SinglC :. Matrix, MultiC :. Vector, MultiC :.  Matrix
+-- SinglC means the input has only one channel, while
+-- MultiC means the input has more than one.
+--
+-- type function composition
+data (f :: * -> *) :. (g :: * -> *) :: * -> *
+-- type function: Identity
+data SinglC :: * -> *
+data MultiC :: * -> *
+
+-- Tags for each form of layer
+data F
+data C
+data A
+data M
+data T (c :: * -> *)
+data S a b
+
+data RunLayer :: * -> * where
+  -- Densely connected layer
+  -- input:   vector of size m
+  -- output:  vector of size n
+  -- weights: matrix of size m x n
+  -- biases:  vector of size n
+  Full :: !(Matrix R) -> !(Vector R) -> RunLayer F
+  -- convolutional layer
+  -- input:  channels of 2D floats, of the same size (a x b), # of input channels:  m
+  -- output: channels of 2D floats, of the same size (c x d), # of output channels: n
+  --         where c = a + 2*padding + 1 - s
+  --               d = b + 2*padding + 1 - t
+  -- feature:  matrix of (s x t), # of features: m x n
+  -- padding:  number of 0s padded at each side of channel
+  -- biases:   bias for each output, # of biases: n
+  Conv  :: !(V.Vector (V.Vector (Matrix R))) -> !(V.Vector R) -> Int -> RunLayer C
+  -- Reshape from channels of matrix to a single vector
+  -- input:  m channels of 2D matrices
+  --         assuming that all matrices are of the same size a x b
+  -- output: 1D vector of the concatenation of all input channels
+  --         its size: m x a x b
+  As1D  :: RunLayer A
+  -- max pooling layer
+  -- input:  channels of 2D floats, of the same size (a x b), # of input channels:  m
+  --         assuming that a and b are both multiple of stride
+  -- output: channels of 2D floats, of the same size (c x d), # of output channels: m
+  --         where c = a / stride
+  --               d = b / stride
+  MaxP :: Int -> RunLayer M
+  -- Activator
+  -- the input can be either a 1D vector, 2D matrix, or channels of either.
+  Activation :: (R->R, R->R) -> RunLayer (T c)
+  -- stacking two components a and b
+  -- the output of a should matches the input of b
+  Stack :: !(RunLayer a) -> !(RunLayer b) -> RunLayer (S a b)
+
+instance Component (RunLayer F) where
+    type Run (RunLayer F) = Identity
+    type Inp (RunLayer F) = Vector R
+    type Out (RunLayer F) = Vector R
+    -- trace is (input, weighted-sum)
+    newtype Trace (RunLayer F) = DTrace (Vector R, Vector R)
+    forwardT (Full w b) !inp =
+        let !bv = (inp <# w) `add` b
+        in return $ DTrace (inp,bv)
+    output (DTrace (_,!a)) = a
+    backward l (DTrace (!iv,!bv)) !odelta rate =
+        let Full w b = l
+            !d = scale (negate rate) odelta
+            !m = iv `outer` d
+            -- back-propagated error at input
+            !idelta = w #> odelta
+            -- update to weights
+            ---- for what reason, could this expression: w `add` (iv `outer` d)
+            ---- entails a huge space leak? especially, neither 'seq' nor
+            ---- 'deepseq' helps a bit. The only workaround is to expand the
+            ---- add function, and call SV.force on the result vector, which
+            ---- explcitly copy and drop reference to orignal computed result.
+            !w'= w `add` m
+            -- !w'= let (r,c) = size w
+            --          dat1 = flatten (tr' w)
+            --          dat2 = flatten (tr' m)
+            --      in matrixFromVector ColumnMajor r c $ SV.force $ dat1 `add` dat2
+            !b'= b `add` d
+            -- !b'= SV.force $ b `add` d
+        in return $ (Full w' b', idelta)
+
+instance Component (RunLayer C) where
+    type Run (RunLayer C) = Identity
+    type Inp (RunLayer C) = V.Vector (Matrix R)
+    type Out (RunLayer C) = V.Vector (Matrix R)
+    -- trace is (input, convoluted output)
+    newtype Trace (RunLayer C) = CTrace (Inp (RunLayer C), V.Vector (Matrix R))
+    forwardT (Conv fs bs p) !inp =
+        let !ov = parallel $ V.zipWith feature
+                               (tr fs) -- feature matrix indexed majorly by each output
+                               bs      -- biases by each output
+        in return $ CTrace (inp,ov)
+      where
+        !osize = let (x,y) = size (V.head inp)
+                     (u,v) = size (V.head $ V.head fs)
+                 in (x+2*p-u+1, y+2*p-v+1)
+        -- transpose the features matrix
+        tr :: V.Vector (V.Vector a) -> V.Vector (V.Vector a)
+        tr uv = let n   = V.length (V.head uv)
+                    !vu = V.map (\i -> V.map (V.! i) uv) $ V.enumFromN 0 n
+                in vu
+        feature :: V.Vector (Matrix R) -> R -> Matrix R
+        feature f b = V.foldl1' add (V.zipWith (layerCorr2 p) f inp) `add` konst b osize
+    output (CTrace (_,a)) = a
+    backward l (CTrace (!iv,!av)) !odelta rate =
+      let Conv fs bs p = l
+          -- update to the feature matrix
+          m :: V.Vector (V.Vector (Matrix R))
+          !m = parallel $ V.zipWith (\flts chn ->
+                            -- chn:  a single input channel
+                            -- flts: all features used for chn
+                            V.zipWith (\f d ->
+                              let upd = scale (negate rate) (layerCorr2 p chn d)
+                              in f `add` upd
+                            ) flts odelta
+                          ) fs iv
+          -- update to the biases
+          b :: V.Vector R
+          !b = V.zipWith (\b d -> b + (negate rate) * sumElements d) bs odelta
+          -- back-propagated error at input
+          idelta :: V.Vector (Matrix R)
+          !idelta = V.map (\f -> V.foldl1' add $ V.zipWith (layerConv2 p) f odelta) fs
+      in --trace ("CL:" ++ show odelta)
+         return $ (Conv m b p, idelta)
+
+instance Component (RunLayer A) where
+  type Run (RunLayer A) = Identity
+  type Inp (RunLayer A) = V.Vector (Matrix R)
+  type Out (RunLayer A) = Vector R
+  -- trace keeps information of (m, axb, b, output)
+  newtype Trace (RunLayer A) = ReshapeTrace (Int, Int, Int, Vector R)
+  forwardT _ !inp =
+    let !b = V.length inp
+        (!r,!c) = size (V.head inp)
+        !o = V.foldr' (\x y -> flatten x SV.++ y) SV.empty inp
+    in return $ ReshapeTrace (b, r*c, c, o)
+  output (ReshapeTrace (_,_,_,a)) = a
+  backward a (ReshapeTrace (b,n,c,_)) !odelta _ =
+    let !idelta = V.fromList $ map (reshape c) $ takesV (replicate b n) odelta
+    in return $ (a, idelta)
+
+instance Component (RunLayer M) where
+  type Run (RunLayer M) = Identity
+  type Inp (RunLayer M) = V.Vector (Matrix R)
+  type Out (RunLayer M) = V.Vector (Matrix R)
+  -- trace is (dimension of pools, index of max in each pool, pooled matrix)
+  -- for each channel.
+  newtype Trace (RunLayer M) = PTrace (V.Vector (IndexOf Matrix, Vector Int, Matrix R))
+  -- forward is to divide the input matrix in stride x stride sub matrices,
+  -- and then find the max element in each sub matrices.
+  forwardT (MaxP stride) !inp = return $ PTrace $ parallel $ V.map mk inp
+    where
+      mk inp = let (!i,!v) = pool stride inp in (size v, i, v)
+  output (PTrace a) = V.map (\(_,_,!o) ->o) a
+  -- use the saved index-of-max in each pool to propagate the error.
+  backward l@(MaxP stride) (PTrace t) odelta _ =
+      let !idelta = V.zipWith gen t odelta in return $ (l, idelta)
+    where
+      gen (!si,!iv,_) od = unpool stride iv od
+
+instance (Component (RunLayer a),
+          Component (RunLayer b),
+          Run (RunLayer a) ~ Identity,
+          Run (RunLayer b) ~ Identity,
+          Out (RunLayer a) ~ Inp (RunLayer b)
+         ) => Component (RunLayer (S a b)) where
+    type Run (RunLayer (S a b)) = Identity
+    type Inp (RunLayer (S a b)) = Inp (RunLayer a)
+    type Out (RunLayer (S a b)) = Out (RunLayer b)
+    newtype Trace (RunLayer (S a b)) = TTrace (Trace (RunLayer b), Trace (RunLayer a))
+    forwardT (Stack a b) !i = do
+        !tra <- forwardT a i
+        !trb <- forwardT b (output tra)
+        return $ TTrace (trb, tra)
+    output (TTrace !a) = output (fst a)
+    backward (Stack a b) (TTrace (!trb,!tra)) !odelta rate = do
+        (b', !odelta') <- backward b trb odelta  rate
+        (a', !idelta ) <- backward a tra odelta' rate
+        return (Stack a' b', idelta)
+
+instance (Container c R) => Component (RunLayer (T (MultiC :. c))) where
+    type Run (RunLayer (T (MultiC :. c))) = Identity
+    type Inp (RunLayer (T (MultiC :. c))) = V.Vector (c R)
+    type Out (RunLayer (T (MultiC :. c))) = V.Vector (c R)
+    newtype Trace (RunLayer (T (MultiC :. c))) = TTraceM (V.Vector (Trace (RunLayer (T (SinglC :. c)))))
+    forwardT (Activation ac) !inp =
+      TTraceM <$> V.mapM (forwardT (Activation ac)) inp
+    output (TTraceM a) = V.map output a
+    backward a@(Activation ac) (TTraceM ts) !odelta r = do
+      idelta <- V.zipWithM (\t d -> snd <$> backward (Activation ac) t d r) ts odelta
+      return (a, idelta)
+
+instance (Container c R) => Component (RunLayer (T (SinglC :. c))) where
+    type Run (RunLayer (T (SinglC :. c))) = Identity
+    type Inp (RunLayer (T (SinglC :. c))) = c R
+    type Out (RunLayer (T (SinglC :. c))) = c R
+    newtype Trace (RunLayer (T (SinglC :. c))) = TTraceS (c R, c R)
+    forwardT (Activation (af,_)) !inp = return $ TTraceS (inp, cmap af inp)
+    output (TTraceS (_,!a)) = a
+    backward a@(Activation (_,ag)) (TTraceS (!iv,_)) !odelta _ = return $ (a, odelta `hadamard` cmap ag iv)
+
+newFLayer :: Int                -- number of input values
+          -> Int                -- number of neurons (output values)
+          -> IO (RunLayer F)    -- new layer
+newFLayer m n =
+    withSystemRandom . asGenIO $ \gen -> do
+        -- we build the weights in column major because in the back-propagation
+        -- algo, the computed update to weights is in column major. So it is
+        -- good for performance to keep the matrix always in column major.
+        w <- buildMatrix (normal 0 0.01 gen) ColumnMajor (m,n)
+        b <- return $ konst 1 n
+        return $ Full w b
+
+newCLayer :: Int                -- number of input channels
+          -> Int                -- number of output channels
+          -> Int                -- size of each feature
+          -> Int                -- size of padding
+          -> IO (RunLayer C)    -- new layer
+newCLayer inpsize outsize sfilter npadding =
+  withSystemRandom . asGenIO $ \gen -> do
+      fs <- V.replicateM inpsize $ V.replicateM outsize $
+              buildMatrix (truncNormal 0 0.1 gen) RowMajor (sfilter, sfilter)
+      bs <- return $ V.replicate outsize 0.1
+      return $ Conv fs bs npadding
+  where
+    truncNormal m s g = do
+      x <- standard g
+      if x >= 2.0 || x <= -2.0
+        then truncNormal m s g
+        else return $! m + s * x
+
+buildMatrix g order (nr, nc) = do
+  vals <- SV.replicateM (nr*nc) (double2Float <$> g)
+  return $ matrixFromVector order nr nc vals
+
+layerCorr2 :: Int -> Matrix R -> Matrix R -> Matrix R
+layerCorr2 p k m = c_corr2d_s k padded
+  where
+    padded = zeroPadded p m
+    (w,_)  = size k
+
+layerConv2 :: Int -> Matrix R -> Matrix R -> Matrix R
+layerConv2 p k m = c_conv2d_s k padded
+  where
+    padded = zeroPadded p m
+    (w,_)  = size k
+
+-- max pool, picking out the maximum element
+-- in each stride x stride sub-matrices.
+-- assuming that the original matrix row and column size are
+-- both multiple of stride
+pool :: Int -> Matrix Float -> (Vector Int, Matrix Float)
+pool 1 mat = let (r,c) = size mat in (SV.replicate (r*c) 0, mat)
+-- pool 2 mat | orderOf mat == RowMajor = c_max_pool2_f mat
+pool stride mat = runST $ do
+  ori <- unsafeThawMatrix mat
+  mxv <- newUndefinedMatrix RowMajor r' c'
+  mxi <- newUndefinedVector (r'*c')
+  forM_ [0..r'-1] $ \i -> do
+    forM_ [0..c'-1] $ \j -> do
+      (n,v) <- unsafeMaxIndEle ori (i*stride) (j*stride) stride stride
+      unsafeWriteVector mxi (i*c'+j) n
+      unsafeWriteMatrix mxv i j v
+  a <- unsafeFreezeVector mxi
+  b <- unsafeFreezeMatrix mxv
+  return (a,b)
+  where
+    (r,c) = size mat
+    r'    = r `div` stride
+    c'    = c `div` stride
+    unsafeMaxIndEle mm x y r c = do
+      mp <- newSTRef 0
+      mv <- newSTRef (-10000.0)
+      forM_ [0..r-1] $ \ i -> do
+        forM_ [0..c-1] $ \ j -> do
+          v1 <- unsafeReadMatrix mm (x+i) (y+j)
+          v0 <- readSTRef mv
+          when (v1 > v0) $ do
+            writeSTRef mv v1
+            writeSTRef mp (i*2+j)
+      p <- readSTRef mp
+      v <- readSTRef mv
+      return (p, v)
+
+-- the reverse of max pool.
+-- assuming idx and mat are of the same size
+unpool :: Int -> Vector Int -> Matrix Float -> Matrix Float
+unpool stride idx mat = runSTMatrix $ do
+  mat' <- newMatrix' 0 r' c'
+  forM_ [0..r-1] $ \i -> do
+    forM_ [0..c-1] $ \j -> do
+      let pos     = idx SV.! (i*c+j)
+      let (oi,oj) = pos `divMod` 2
+      let val     = mat `atIndex` (i,j)
+      unsafeWriteMatrix mat' (i*stride+oi) (j*stride+oj) val
+  return mat'
+  where
+    (r,c) = size mat
+    (r',c') = (r*stride, c*stride)
+
+-- a slightly faster way to pading the matrix
+-- camparing to fromBlocks provided by hmatrix.
+zeroPadded :: Int -> Matrix Float -> Matrix Float
+zeroPadded p mat = runSTMatrix $ do
+  mat' <- newMatrix' 0 r' c'
+  setMatrix mat' p p mat
+  return mat'
+  where
+    (r,c) = size mat
+    (r',c') = (r+2*p,c+2*p)
+
+-- a slightly faster version of newMatrix, which based
+-- directly on lower level Vector.Storable creation.
+newMatrix' :: SVM.Storable t => t -> Int -> Int -> ST s (STMatrix s t)
+newMatrix' v r c = do
+  vec <- SVM.replicate (r*c) v
+  vec <- SV.unsafeFreeze vec
+  unsafeThawMatrix $ reshape c vec
diff --git a/Data/NeuralNetwork/Backend/HMatrix/Utils.hs b/Data/NeuralNetwork/Backend/HMatrix/Utils.hs
new file mode 100644
--- /dev/null
+++ b/Data/NeuralNetwork/Backend/HMatrix/Utils.hs
@@ -0,0 +1,177 @@
+{-# LANGUAGE BangPatterns, FlexibleInstances, FlexibleContexts, ForeignFunctionInterface #-}
+module Data.NeuralNetwork.Backend.HMatrix.Utils where
+import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.Devel
+import Numeric.GSL.Fourier
+import Data.Complex
+import Control.Exception
+import Control.DeepSeq
+import Control.Parallel
+import Control.Parallel.Strategies
+import Control.Monad
+import qualified Data.Vector as VecB
+import qualified Data.Vector.Generic as VecGeneric
+import qualified Data.Vector.Fusion.Bundle as VecFusion
+import qualified Data.Vector.Fusion.Bundle.Monadic as VecFusionM
+import qualified Data.Vector.Storable as VecS
+import qualified Data.Vector.Storable.Mutable as VecM
+import Foreign.Ptr ( Ptr )
+import Foreign.C.Types ( CInt(..) )
+import System.IO.Unsafe ( unsafePerformIO )
+
+-- fft2d :: Matrix (Complex Double) -> Matrix (Complex Double)
+-- fft2d m = let !x = fromRows $ map fft $ unsafeToRows m
+--               !y = fromColumns $ map fft $ toColumns x
+--           in y
+-- ifft2d :: Matrix (Complex Double) -> Matrix (Complex Double)
+-- ifft2d m = let !x = fromRows $ map ifft $ unsafeToRows m
+--                !y = fromColumns $ map ifft $ toColumns x
+--            in y
+
+-- fft2d :: Matrix (Complex Double) -> Matrix (Complex Double)
+-- fft2d m = let rh:rr = map fft $ force $ toRows m
+--               x = fromRows $ withStrategy (parList rdeepseq) rr `pseq` (rh : rr)
+--               sh:ss = map fft $ force $ toColumns x
+--               y = fromColumns $ withStrategy (parList rdeepseq) ss `pseq` (sh : ss)
+--           in y
+--
+-- ifft2d :: Matrix (Complex Double) -> Matrix (Complex Double)
+-- ifft2d m = let rh:rr = map fft $ force $ toRows m
+--                x = fromRows $ withStrategy (parList rdeepseq) rr `pseq` (rh : rr)
+--                sh:ss = map fft $ force $ toColumns x
+--                y = fromColumns $ withStrategy (parList rdeepseq) ss `pseq` (sh : ss)
+--            in y
+
+-- conv2d_b :: (Numeric t, ConvFD t, Container Vector t, Container Matrix t)
+--          => Matrix t -> Matrix t -> Matrix t
+-- conv2d_b !k !m | w1 > w2 && h1 < h2 = error "convolution cannot be performed"
+--                | w1 < w2 && h1 > h2 = error "convolution cannot be performed"
+--                | w1 > w2   = conv2d_b m k
+--                | otherwise =
+--     -- convolution via FFT is actually cyclic conv. so we need to 0-pad the
+--     -- matrix being convoluted by half the size of the kernel matrix.
+--     -- the kernel matrix is also 0-padded to be the equal size.
+--     -- finall, extra
+--     let !hw = w1 `div` 2
+--         !hh = h1 `div` 2
+--         !z1 = konst 0 (w2-w1+hw,h2-h1+hh)
+--         !z2 = konst 0 (hw, hh)
+--         m1' = fft2d $ fromBlocks [[m1,0],[0,z1]]
+--         m2' = fft2d $ fromBlocks [[m2,0],[0,z2]]
+--         mr  = ifft2d $ (force m1' `par` (force m2' `pseq` hadamard m1' m2'))
+--         ms  = subMatrix (w1-1,h1-1) (w2-w1+1,h2-h1+1) $ fst . fromComplex $ mr
+--     in fromDouble $ ms
+--   where
+--      !m1 = complex $ toDouble k
+--      !m2 = complex $ toDouble m
+--      (w1,h1) = size m1
+--      (w2,h2) = size m2
+--
+-- corr2d_b k m | w > s && h < t = error "convolution cannot be performed"
+--              | w < s && h > t = error "convolution cannot be performed"
+--              | w > s     = conv2d_b (rotate m) k
+--              | otherwise = conv2d_b (rotate k) m
+--   where
+--     (w,h)   = size k
+--     (s,t)   = size m
+--
+-- corr2d_s :: (Numeric t, Container Vector t, Container Matrix t)
+--          => Matrix t -> Matrix t -> Matrix t
+-- corr2d_s k m | w > s && h < t = error "correlation cannot be performed"
+--              | w < s && h > t = error "correlation cannot be performed"
+--              | w > s = corr2d_s m k
+--              | otherwise = final
+--   where
+--     (w,h) = size k
+--     (s,t) = size m
+--     (u,v) = (s-w+1, t-h+1)
+--     subs  = map (\s->subMatrix s (w,h) m) $ [(x,y) | x<-[0..u-1], y<-[0..v-1]]
+--     -- use unsafe* methods to create the intermediate matrix fast.
+--     t_rows = u*v
+--     t_cols = w*h
+--     !transformed = matrixFromVector RowMajor t_rows t_cols $ VecS.create $ do
+--         mat <- VecM.new (t_rows*t_cols)
+--         forM_ (zip [0..] subs) $ \(ri,rm) -> do
+--             let bs = t_cols*ri
+--             forM_ (zip [0..] $ unsafeToRows rm) $ \(ci, rv) -> do
+--                 let tv = VecM.unsafeSlice (bs+h*ci) h mat
+--                 {-# SCC "corr-data-copy" #-} VecS.unsafeCopy tv rv
+--         return mat
+--     !weights = flatten k
+--     !final   = reshape v $ transformed #> weights
+--
+-- conv2d_s k m | w > s && h < t = error "convolution cannot be performed"
+--              | w < s && h > t = error "convolution cannot be performed"
+--              | w > s     = corr2d_s (rotate m) k
+--              | otherwise = corr2d_s (rotate k) m
+--   where
+--     (w,h)   = size k
+--     (s,t)   = size m
+
+-- class ConvFD t where
+--   fromDouble :: Matrix Double -> Matrix t
+--   toDouble   :: Matrix t -> Matrix Double
+-- instance ConvFD Double where
+--   fromDouble = id
+--   toDouble = id
+-- instance ConvFD Float where
+--   fromDouble = single
+--   toDouble = double
+--
+-- rotate :: Element t => Matrix t -> Matrix t
+-- rotate = fliprl . flipud
+
+foreign import ccall unsafe corr_sf_general ::
+  CInt ->
+  CInt -> CInt -> CInt -> Ptr Float ->
+  CInt -> CInt -> CInt -> Ptr Float ->
+  Ptr Float -> IO CInt
+
+data CConvType = CConv | CCorr
+
+c_corr2d_g :: CConvType -> Matrix Float -> Matrix Float -> Matrix Float
+c_corr2d_g y k m | w > s = c_corr2d_s m k
+                 | orderOf k == RowMajor && orderOf m == RowMajor =
+                   let (r,c) = (s-w+1,t-h+1)
+                       v = unsafePerformIO $ do
+                             v <- VecM.unsafeNew (r * c)
+                             VecM.unsafeWith v $ \rp ->
+                               apply k id $ \kr kc ks kt kp ->
+                                 apply m id $ \mr mc ms mt mp ->
+                                   case y of
+                                     CCorr -> corr_sf_general 0 kr kc ks kp mr mc ms mp rp
+                                     CConv -> corr_sf_general 1 kr kc ks kp mr mc ms mp rp
+                             VecS.unsafeFreeze v
+                   in matrixFromVector RowMajor r c v
+                | otherwise = error "column major matrix not supported"
+  where
+    (w,h) = size k
+    (s,t) = size m
+
+c_corr2d_s = c_corr2d_g CCorr
+c_conv2d_s = c_corr2d_g CConv
+
+-- parallel !vec = vec
+parallel :: NFData a => VecB.Vector a -> VecB.Vector a
+parallel vec = (VecB.tail vec `using` parvec) `pseq` vec
+  where
+    parvec = VecB.mapM (rparWith rdeepseq)
+
+-- foreign import ccall unsafe pool2_f :: CInt -> CInt -> CInt -> Ptr Float ->
+--                                        Ptr Float -> Ptr CInt -> IO ()
+-- c_max_pool2_f :: Matrix Float -> (Vector Int, Matrix Float)
+-- c_max_pool2_f mat
+--   | orderOf mat == RowMajor = unsafePerformIO $ do
+--         ind <- VecM.unsafeNew (r' * c')
+--         mx  <- VecM.unsafeNew (r' * c')
+--         VecM.unsafeWith ind $ \pind ->
+--           VecM.unsafeWith mx $ \pmax ->
+--             apply mat id $ \mr mc ms mt mp ->
+--               pool2_f mr mc ms mp pmax pind
+--         ind <- VecS.unsafeFreeze ind
+--         mx  <- VecS.unsafeFreeze mx
+--         return (VecS.map fromIntegral ind, matrixFromVector RowMajor r' c' mx)
+--   where
+--     (r,c) = size mat
+--     r'    = r `div` 2
+--     c'    = c `div` 2
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,29 @@
+BSD 3-Clause License
+
+Copyright (c) 2016, Jiasen Wu
+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.
+
+* Neither the name of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+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 HOLDER 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/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/cbits/conv.c b/cbits/conv.c
new file mode 100644
--- /dev/null
+++ b/cbits/conv.c
@@ -0,0 +1,70 @@
+#include <stdlib.h>
+#include <string.h>
+#include "cblas.h"
+
+#define AT(p,s,x,y) ((p)+(s)*(x)+(y))
+
+
+// 2M byte working storage per thread.
+// sufficient to handle correlation/convolution wthin
+// the size (5x5, 128x128) or (7x7, 100x100)
+#define WORKINGSTORAGE 524288
+static _Thread_local float workingSto[WORKINGSTORAGE];
+
+/*
+input:  two matrices, assuming in row major, may not be in continuous. i.e.
+        one row has only column number of elements, but its length is of
+        the stride.
+        mat1 is the kernel
+        mat2 is the source
+        assuming that row1<=row2 and col1 <= col2
+output: continuous row-major matrix of size u x v
+        where u = row2-row1+1
+              v = col2-co11+1
+*/
+
+int corr_sf_general(
+    int reversemat1, int row1, int col1, int stride1, float *mat1,
+    int row2, int col2, int stride2, float *mat2, float *mat3)
+{
+    int u = row2-row1+1;
+    int v = col2-col1+1;
+    if ((u*v+1)*row1*col1 > WORKINGSTORAGE) return -1;
+    float *ws = workingSto;
+    for(int i=0;i<u;i++) {
+        for(int j=0;j<v;j++) {
+            for(int k=0;k<row1;k++) {
+                float *src=AT(mat2,stride2,i+k,j);
+                memcpy(ws, src, col1*sizeof(float));
+                ws += col1;
+            }
+        }
+    }
+    float *vw = mat1;
+    if (stride1 > col1 || reversemat1) {
+        // we have expecting an extra row1*col1 elements at hand, pointed
+        // by ws at this point of time, and used for a continues
+        // copy of the kernel matrix.
+        vw = ws;
+        float *p1=vw, *p2=mat1;
+        for(int i=0;i<row1;i++) {
+            memcpy(p1,p2,col1*sizeof(float));
+            p1 += col1;
+            p2 += stride1;
+        }
+    }
+    // reverse the kernel when do convolution.
+    if(reversemat1) {
+        int sz = row1*col1;
+        for(int i=0;i<sz/2;i++) {
+            float t = vw[i];
+            vw[i] = vw[sz-1-i];
+            vw[sz-1-i] = t;
+        }
+    }
+    // execute the matrix vector multiplication
+    ws = workingSto;
+    cblas_sgemv(CblasRowMajor, CblasNoTrans, u*v, row1*col1,
+    		    1.0, ws, row1*col1, vw, 1, 0, mat3, 1);
+    return 0;
+}
diff --git a/neural-network-hmatrix.cabal b/neural-network-hmatrix.cabal
new file mode 100644
--- /dev/null
+++ b/neural-network-hmatrix.cabal
@@ -0,0 +1,25 @@
+name:                neural-network-hmatrix
+version:             0.1.0.0
+license-file:        LICENSE
+license:             BSD3
+author:              Jiasen Wu
+maintainer:          jiasenwu@hotmail.com
+Category:            AI
+Synopsis:            Yet Another High Performance and Extendable Neural Network in Haskell
+Description:         Provides execution backend of neural network on top of hmatrix.
+build-type:          Simple
+cabal-version:       >=1.10
+
+library
+  exposed-modules:     Data.NeuralNetwork.Backend.HMatrix
+  other-modules:       Data.NeuralNetwork.Backend.HMatrix.Utils
+                       Data.NeuralNetwork.Backend.HMatrix.Layers
+  build-depends:       base >= 4.7 && < 5, hmatrix, hmatrix-gsl, mwc-random, mtl,
+                       vector, deepseq, parallel, neural-network-base
+  default-language:    Haskell2010
+  C-sources:           cbits/conv.c
+  cc-options:          -O3 -std=c11 -march=native
+  if os(windows)
+      extra-libraries: openblas
+  if os(linux)
+      extra-libraries: blas