packages feed

neural-network-blashs (empty) → 0.1.0.0

raw patch · 13 files changed

+1700/−0 lines, 13 filesdep +QuickCheckdep +basedep +blas-hsbuild-type:Customsetup-changed

Dependencies added: QuickCheck, base, blas-hs, constraints, ghc-prim, hmatrix, hspec, mtl, mwc-random, neural-network-base, vector

Files

+ Data/NeuralNetwork/Backend/BLASHS.hs view
@@ -0,0 +1,117 @@+------------------------------------------------------------
+-- |
+-- Module      :  Data.NeuralNetwork.Backend.BLASHS
+-- Description :  A backend for neural network on top of 'blas-hs'
+-- Copyright   :  (c) 2016 Jiasen Wu
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Jiasen Wu <jiasenwu@hotmail.com>
+-- Stability   :  experimental
+-- Portability :  portable
+--
+--
+-- This module supplies a backend for the neural-network-base
+-- package. This backend is implemented on top of the blas-hs
+-- package and optimised with SIMD.
+------------------------------------------------------------
+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeFamilies #-}
+module Data.NeuralNetwork.Backend.BLASHS (
+  -- module Data.NeuralNetwork.Backend.BLASHS.Layers,
+  module Data.NeuralNetwork.Backend.BLASHS.Utils,
+  ByBLASHS(..),
+  ErrCode(..),
+  cost'
+) where
+
+import Data.NeuralNetwork hiding (relu, relu', cost')
+import Data.NeuralNetwork.Backend.BLASHS.Layers
+import Data.NeuralNetwork.Backend.BLASHS.Utils
+import Data.NeuralNetwork.Backend.BLASHS.SIMD
+import Control.Monad.Except
+import Data.Constraint (Dict(..))
+
+-- | Compilation of the specification of a neural network is carried out in
+-- the 'Err' monad, and the possible errors are characterized by 'ErrCode'.
+type Err     = ExceptT ErrCode IO
+data ErrCode = ErrMismatch
+
+-- | The backend data type
+data ByBLASHS = ByBLASHS
+
+-- | Neural network specified to start with 1D / 2D input
+instance (HeadSize z, TranslateBody s,
+          Component (RunLayer (SpecToTag s)),
+          Run (RunLayer (SpecToTag s)) ~ IO)
+       => Backend ByBLASHS (z :++ s) where
+  type Env ByBLASHS = Err
+  type ConvertFromSpec (z :++ s) = RunLayer (SpecToTag s)
+  compile _ (a :++ l)= trans (hsize a) l
+  witness _ _ = Dict
+
+instance RunInEnv IO Err where
+  run = liftIO
+
+-- It is necessary to propagate the size along the layers,
+-- because fullconnect and convolution need to know
+-- the previous size.
+data LayerSize = D1 Int | D2 Int Int Int
+
+-- 'HeadSize' is class for the input layer
+class HeadSize l where
+  hsize :: l -> LayerSize
+instance HeadSize SpecIn1D where
+  hsize (In1D n) = D1 n
+instance HeadSize SpecIn2D where
+  hsize (In2D m n) = D2 1 m n
+-- 'BodySize' is class for the actual computational layers
+class BodySize l where
+  bsize :: LayerSize -> l -> LayerSize
+instance BodySize SpecReshape2DAs1D where
+  bsize (D2 k m n) _ = D1 (k*m*n)
+instance BodySize SpecFullConnect where
+  bsize _ (FullConnect n) = D1 n
+instance BodySize SpecConvolution where
+  bsize (D2 _ m n) (Convolution k f p) = D2 k (m+2*p-f+1) (n+2*p-f+1)
+instance BodySize SpecMaxPooling where
+  bsize (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 :: LayerSize -> s -> Err (RunLayer (SpecToTag s))
+
+instance TranslateBody SpecFullConnect where
+  -- 'SpecFullConnect' is translated to a two-layer component
+  -- a full-connect, followed by a relu activation (1D, single channel)
+  type SpecToTag SpecFullConnect = S F (T SinglVec)
+  trans (D1 s) (FullConnect n) = do u <- lift $ newFLayer s n
+                                    return $ Stack u (Activation (relu, relu'))
+  trans _ _ = throwError ErrMismatch
+
+instance TranslateBody SpecConvolution where
+  -- 'SpecConvolution' is translated to a two-layer component
+  -- a convolution, following by a relu activation (2D, multiple channels)
+  type SpecToTag SpecConvolution = S C (T MultiMat)
+  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 SpecMaxPooling where
+  -- 'MaxPooling' is translated to a max-pooling component.
+  type SpecToTag SpecMaxPooling = P
+  trans (D2 _ _ _) (MaxPooling n) = return (MaxP n)
+  trans (D1 _)     _              = throwError ErrMismatch
+
+instance TranslateBody SpecReshape2DAs1D where
+  -- 'SpecReshape2DAs1D' is translated to a reshaping component.
+  type SpecToTag SpecReshape2DAs1D = A
+  trans (D2 _ _ _) _ = return As1D
+  trans (D1 _)     _ = throwError ErrMismatch
+
+instance (TranslateBody a, TranslateBody c, BodySize a) => TranslateBody (a :++ c) where
+  -- ':++' is translated to the stacking component.
+  type SpecToTag (a :++ b) = S (SpecToTag a) (SpecToTag b)
+  trans s (a :++ c) = do u <- trans s a
+                         v <- trans (bsize s a) c
+                         return $ Stack u v
+ Data/NeuralNetwork/Backend/BLASHS/Layers.hs view
@@ -0,0 +1,285 @@+------------------------------------------------------------
+-- |
+-- Module      :  Data.NeuralNetwork.Backend.BLASHS.Utils
+-- Description :  A backend for neuralnetwork with blas-hs.
+-- Copyright   :  (c) 2016 Jiasen Wu
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Jiasen Wu <jiasenwu@hotmail.com>
+-- Stability   :  experimental
+-- Portability :  portable
+--
+--
+-- This module supplies a high level abstraction of the rather
+-- low-level blas-hs interfaces.
+------------------------------------------------------------
+{-# LANGUAGE BangPatterns, TypeFamilies, TypeOperators, FlexibleInstances, FlexibleContexts, GADTs #-}
+module Data.NeuralNetwork.Backend.BLASHS.Layers(
+  SinglVec, MultiMat, F, C, A, P, T, S, RunLayer(..),
+  newFLayer, newCLayer
+) where
+
+import qualified Data.Vector as V
+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.NeuralNetwork
+import Data.NeuralNetwork.Backend.BLASHS.Utils
+import Data.NeuralNetwork.Backend.BLASHS.SIMD
+
+type R = Float
+type M = IO
+
+-- | We parameterise the activation layer T, where the parameter indicates how
+-- elements are contained:
+data SinglVec
+data MultiMat
+
+-- | tag for the full-connect component
+data F
+-- | tag for the convolution component
+data C
+-- | tag for the component that converts 2D as 1D
+data A
+-- | tag for the max-pooling component
+data P
+-- | tag for the activation component
+data T c
+-- | tag for the stacking component
+data S a b
+
+-- | basic components of neural network
+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 :: !(DenseMatrix R) -> !(DenseVector 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 (DenseMatrixArray 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 P
+  -- | Activator
+  -- the input can be either a 1D vector, 2D matrix, or channels of either.
+  Activation :: (SIMDPACK R -> SIMDPACK R, SIMDPACK R -> SIMDPACK 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) = IO
+    type Inp (RunLayer F) = DenseVector R
+    type Out (RunLayer F) = DenseVector R
+    -- trace is (input, weighted-sum)
+    newtype Trace (RunLayer F) = DTrace (DenseVector R, DenseVector R)
+    forwardT (Full !w !b) !inp = do
+        bv <- newDenseVectorCopy b
+        bv <<+ inp :<# w
+        return $ DTrace (inp,bv)
+    output (DTrace (_,!a)) = a
+    backward (Full !w !b) (DTrace (!iv,!bv)) !odelta rate = do
+        -- back-propagated error at input
+        idelta <- newDenseVector (fst $ size w)
+        idelta <<= w :#> odelta
+        -- odelta is not used any more, so we reuse it for an intermediate value.
+        odelta <<= Scale (negate rate)
+        w <<+ iv :## odelta
+        b <<= b  :.+ odelta
+        return (Full w b, idelta)
+
+instance Component (RunLayer A) where
+  type Run (RunLayer A) = IO
+  type Inp (RunLayer A) = V.Vector (DenseMatrix R)
+  type Out (RunLayer A) = DenseVector R
+  -- trace keeps information of (m, a, b, output)
+  newtype Trace (RunLayer A) = ReshapeTrace (Int, Int, Int, DenseVector R)
+  forwardT _ !inp = do
+    let !b = V.length inp
+        (!r,!c) = size (V.head inp)
+    o <- denseVectorConcat $ V.map m2v inp
+    return $ ReshapeTrace (b, r, c, o)
+  output (ReshapeTrace (_,_,_,a)) = a
+  backward a (ReshapeTrace (b,r,c,_)) !odelta _ = do
+    let !idelta = V.map (v2m r c) $ denseVectorSplit b (r*c) odelta
+    return $ (a, idelta)
+
+instance Component (RunLayer C) where
+  type Run (RunLayer C) = IO
+  type Inp (RunLayer C) = V.Vector (DenseMatrix R)
+  type Out (RunLayer C) = V.Vector (DenseMatrix R)
+  -- trace is (input, convoluted output)
+  newtype Trace (RunLayer C) = CTrace (Inp (RunLayer C), Out (RunLayer C))
+  forwardT (Conv fss bs pd) !inp = do
+    ma <- newDenseMatrixArray outn outr outc
+    V.zipWithM_ (\fs i -> corr2 pd (denseMatrixArrayToVector fs) i (ma <<+)) fss inp
+    let ov = denseMatrixArrayToVector ma
+    V.zipWithM_ (\m b -> m <<= Apply (plus (konst b))) ov bs
+    return $ CTrace (inp, ov)
+    where
+      outn = V.length bs
+      (outr,outc) = let (x,y)   = size (V.head inp)
+                        (_,u,v) = size (V.head fss)
+                    in (x+2*pd-u+1, y+2*pd-v+1)
+  output (CTrace (_,!a)) = a
+  backward (Conv fss bs pd) (CTrace (iv, av)) !odelta rate = do
+    let (ir,ic) = size (V.head iv)
+    idelta <- newDenseMatrixArray (V.length iv) ir ic
+    fss'   <- transpose fss
+    -- a + 2p - k + 1 = b
+    -- b + 2q - a + 1 = a
+    -- -------------------
+    --    q = k - p
+    -- where
+    --   a = |i|, k = |f|, b = |o|
+    let qd = let (kr,_) = size (V.head $ V.head fss') in kr-1-pd
+    V.zipWithM_ (\fs d -> conv2 qd fs d (idelta <<+)) fss' odelta
+    !nb <- V.zipWithM (\b d -> do s <- sumElements d
+                                  return $ b + negate rate * s
+                      ) bs odelta
+    -- when updating kernels, it originally should be
+    -- conv2 pd iv od. But we use the equalivalent form
+    -- conv2 (|od|-|iv|+pd) od iv. Because there are typically
+    -- more output channels than input.
+    -- a + 2p - b + 1 = c
+    -- b + 2q - a + 1 = c
+    -- ------------------
+    --    q = a - b + p
+    -- where
+    --  a = |o|, b = |i|, c = |f|
+    let qd = let (or,_) = size (V.head odelta) in or - ir + pd
+    V.zipWithM_ (\fs i -> do
+                  -- i:  one input channel
+                  -- fs: all features used for chn
+                  corr2 qd odelta i ((fs <<+) . Scale' (negate rate))
+                ) fss iv
+    let !ideltaV = denseMatrixArrayToVector idelta
+    return $ (Conv fss nb pd, ideltaV)
+instance (Component (RunLayer a),
+          Component (RunLayer b),
+          Run (RunLayer a) ~ IO,
+          Run (RunLayer b) ~ IO,
+          Out (RunLayer a) ~ Inp (RunLayer b)
+         ) => Component (RunLayer (S a b)) where
+    type Run (RunLayer (S a b)) = IO
+    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)) !odeltb rate = do
+        (b', !odelta) <- backward b trb odeltb rate
+        (a', !idelta) <- backward a tra odelta rate
+        return (Stack a' b', idelta)
+
+instance Component (RunLayer (T SinglVec)) where
+    type Run (RunLayer (T SinglVec)) = IO
+    type Inp (RunLayer (T SinglVec)) = DenseVector R
+    type Out (RunLayer (T SinglVec)) = DenseVector R
+    newtype Trace (RunLayer (T SinglVec)) = TTraceS (DenseVector R, DenseVector R)
+    forwardT (Activation (af,_)) !inp = do
+      out <- newDenseVectorCopy inp
+      out <<= Apply af
+      return $ TTraceS (inp, out)
+    output (TTraceS (_,!a)) = a
+    backward a@(Activation (_,ag)) (TTraceS (!iv,_)) !odelta _ = do
+      idelta <- newDenseVectorCopy iv
+      idelta <<= Apply ag
+      idelta <<= odelta :.* idelta
+      return $ (a, idelta)
+
+instance Component (RunLayer (T MultiMat)) where
+    type Run (RunLayer (T MultiMat)) = IO
+    type Inp (RunLayer (T MultiMat)) = V.Vector (DenseMatrix R)
+    type Out (RunLayer (T MultiMat)) = V.Vector (DenseMatrix R)
+    newtype Trace (RunLayer (T MultiMat)) = TTraceM (V.Vector (DenseMatrix R), V.Vector (DenseMatrix R))
+    forwardT (Activation (af,_)) !inp = do
+      out <- V.mapM (\i -> do o <- newDenseMatrixCopy i
+                              o <<= Apply af
+                              return o
+                    ) inp
+      return $ TTraceM (inp, out)
+    output (TTraceM (_,!a)) = a
+    backward a@(Activation (_,ag)) (TTraceM (!iv,_)) !odelta _ = do
+      idelta <- V.zipWithM (\i d -> do o <- newDenseMatrixCopy i
+                                       o <<= Apply ag
+                                       o <<= d :.* o
+                                       return o
+                           ) iv odelta
+      return $ (a, idelta)
+
+instance Component (RunLayer P) where
+  type Run (RunLayer P) = IO
+  type Inp (RunLayer P) = V.Vector (DenseMatrix R)
+  type Out (RunLayer P) = V.Vector (DenseMatrix R)
+  -- trace is (dimension of pools, index of max in each pool, pooled matrix)
+  -- for each channel.
+  newtype Trace (RunLayer P) = PTrace (V.Vector ((Int,Int), DenseVector Int, DenseMatrix 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 = V.mapM mk inp >>= return . PTrace
+    where
+      mk inp = do
+        (!i,!v) <- pool stride inp
+        return (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 _ = do
+      !idelta <- V.zipWithM gen t odelta
+      return $ (l, idelta)
+    where
+      gen (!si,!iv,_) od = unpool stride iv od
+
+-- | create a new full connect component
+newFLayer :: Int                -- ^ number of input values
+          -> Int                -- ^ number of neurons (output values)
+          -> IO (RunLayer F)    -- ^ the new layer
+newFLayer m n =
+    withSystemRandom . asGenIO $ \gen -> do
+        raw <- newDenseVectorByGen (double2Float <$> normal 0 0.01 gen) (m*n)
+        let w = v2m m n raw
+        b <- newDenseVectorConst n 1
+        return $ Full w b
+
+-- | create a new convolutional component
+newCLayer :: Int                -- ^ number of input channels
+          -> Int                -- ^ number of output channels
+          -> Int                -- ^ size of each feature
+          -> Int                -- ^ size of padding
+          -> IO (RunLayer C)    -- ^ the new layer
+newCLayer inpsize outsize sfilter npadding =
+  withSystemRandom . asGenIO $ \gen -> do
+      fss <- V.replicateM inpsize $ do
+              raw <- newDenseVectorByGen (double2Float <$> truncNormal 0 0.1 gen) (outsize*sfilter*sfilter)
+              return $ v2ma outsize sfilter sfilter raw
+      bs <- return $ V.replicate outsize 0.1
+      return $ Conv fss 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
+ Data/NeuralNetwork/Backend/BLASHS/Utils.hs view
@@ -0,0 +1,518 @@+------------------------------------------------------------
+-- |
+-- Module      :  Data.NeuralNetwork.Backend.BLASHS.Utils
+-- Description :  A backend for neuralnetwork with blas-hs.
+-- Copyright   :  (c) 2016 Jiasen Wu
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Jiasen Wu <jiasenwu@hotmail.com>
+-- Stability   :  experimental
+-- Portability :  portable
+--
+--
+-- This module supplies a high level abstraction of the rather
+-- low-level blas-hs interfaces.
+------------------------------------------------------------
+{-# LANGUAGE TypeFamilies, TypeOperators, GADTs #-}
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+{-# LANGUAGE BangPatterns #-}
+module Data.NeuralNetwork.Backend.BLASHS.Utils (
+  DenseVector(..),
+  DenseMatrix(..),
+  DenseMatrixArray(..),
+  newDenseVector,
+  newDenseVectorCopy,
+  newDenseVectorConst,
+  newDenseVectorByGen,
+  newDenseMatrix,
+  newDenseMatrixConst,
+  newDenseMatrixCopy,
+  newDenseMatrixArray,
+  Size(..),
+  denseVectorToVector,
+  denseVectorConcat,
+  denseVectorSplit,
+  denseMatrixArrayAt,
+  denseMatrixArrayToVector,
+  denseMatrixArrayFromVector,
+  v2m, m2v, v2ma, ma2v,
+  Op(..), AssignTo(..),
+  sumElements, corr2, conv2, pool, unpool, transpose
+) where
+
+import Blas.Generic.Unsafe
+import Blas.Primitive.Types
+import qualified Data.Vector as BV
+import qualified Data.Vector.Storable as SV
+import qualified Data.Vector.Storable.Mutable  as V
+import qualified Data.Vector.Storable.Internal as V
+import Control.Exception
+import Control.Monad
+import Data.IORef
+import Foreign.Marshal.Array (advancePtr)
+import Data.NeuralNetwork.Backend.BLASHS.SIMD
+
+-- | mutable vector type
+newtype DenseVector a = DenseVector (V.IOVector a)
+
+-- | mutable matrix type (row-major)
+data DenseMatrix a = DenseMatrix {-# UNPACK #-}!Int {-# UNPACK #-}!Int {-# UNPACK #-}!(V.IOVector a)
+
+-- | array of DenseMatrix, which are identical in size.
+data DenseMatrixArray a = DenseMatrixArray {-# UNPACK #-}!Int {-# UNPACK #-}!Int {-# UNPACK #-}!Int {-# UNPACK #-}!(V.IOVector a)
+
+-- | create a new 'DenseVector'
+newDenseVector :: V.Storable a => Int -> IO (DenseVector a)
+newDenseVector sz = DenseVector <$> V.new sz
+
+-- | create a copy 'DenseVector' from another
+newDenseVectorCopy :: V.Storable a => DenseVector a -> IO (DenseVector a)
+newDenseVectorCopy (DenseVector v) = V.clone v >>= return . DenseVector
+
+-- | create a new 'DenseVector' of some constant
+newDenseVectorConst:: V.Storable a => Int -> a -> IO (DenseVector a)
+newDenseVectorConst n v = V.replicate n v >>= return . DenseVector
+
+-- | create a new 'DenseVector' by a random generator
+newDenseVectorByGen :: V.Storable a => IO a -> Int -> IO (DenseVector a)
+newDenseVectorByGen g n = do
+  vals <- V.replicateM n g
+  return $ DenseVector vals
+
+-- | create a new 'DenseMatrix'
+newDenseMatrix :: V.Storable a => Int -- ^ number of rows
+                               -> Int -- ^ number of columns
+                               -> IO (DenseMatrix a)
+newDenseMatrix r c = DenseMatrix r c <$> V.new (r*c)
+
+-- | create a new 'DenseMatrix' of some constant
+newDenseMatrixConst:: V.Storable a => Int -> Int -> a -> IO (DenseMatrix a)
+newDenseMatrixConst r c v = V.replicate (r*c) v >>= return . DenseMatrix r c
+
+-- | create a copy 'DenseMatrix' from another
+newDenseMatrixCopy :: V.Storable a => DenseMatrix a -> IO (DenseMatrix a)
+newDenseMatrixCopy (DenseMatrix r c v) = V.clone v >>= return . DenseMatrix r c
+
+-- | create a new 'DenseMatrixArray'
+newDenseMatrixArray :: V.Storable a => Int -- ^ number of DenseMatrix
+                                    -> Int -- ^ number of rows
+                                    -> Int -- ^ number of columns
+                                    -> IO (DenseMatrixArray a)
+newDenseMatrixArray n r c = DenseMatrixArray n r c <$> V.new (n*r*c)
+
+-- | get the 'DenseMatrix' from 'DenseMatrixArray' at some position
+denseMatrixArrayAt :: V.Storable a => DenseMatrixArray a -> Int -> DenseMatrix a
+denseMatrixArrayAt (DenseMatrixArray n r c v) i =
+  assert (i >= 0 && i < n) $ let seg = r*c in DenseMatrix r c (V.unsafeSlice (i*seg) seg v)
+
+-- | convert 'DenseMatrixArray' to a vector of 'DenseMatrix' (no copy)
+denseMatrixArrayToVector :: V.Storable a => DenseMatrixArray a -> BV.Vector (DenseMatrix a)
+denseMatrixArrayToVector (DenseMatrixArray n r c v) =
+  let seg = r*c in BV.fromList [DenseMatrix r c (V.unsafeSlice (i*seg) seg v) | i <- [0..n-1]]
+
+-- | convert a vector of 'DenseMatrix' to 'DenseMatrixArray'
+-- If all the matrices are orignally placed consecutively in storage, the result
+-- is simply a type-cast. Otherwise, a new storage is obtained, and matrices are
+-- copied.
+denseMatrixArrayFromVector :: V.Storable a => BV.Vector (DenseMatrix a) -> IO (DenseMatrixArray a)
+denseMatrixArrayFromVector vm = do
+  let n = BV.length vm
+      DenseMatrix r c (V.MVector _ ptr0) = BV.head vm
+  DenseVector raw <- denseVectorConcat (BV.map m2v vm)
+  return $ DenseMatrixArray n r c raw
+
+-- | type cast from 'DenseVector' to 'DenseMatrix'
+v2m r c (DenseVector v) = DenseMatrix r c v
+-- | type cast from 'DenseMatrix' to 'DenseVector'
+m2v (DenseMatrix _ _ v) = DenseVector v
+-- | type cast from 'DenseVector' to 'DenseMatrixArray'
+v2ma n r c (DenseVector v) = assert (V.length v == n*r*c) $ DenseMatrixArray n r c v
+-- | type cast from 'DenseMatrixArray' to 'DenseVector'
+ma2v (DenseMatrixArray n r c v) = DenseVector v
+
+-- | convert a 'DenseVector' to a vector of elements
+denseVectorToVector :: V.Storable a => DenseVector a -> IO (BV.Vector a)
+denseVectorToVector (DenseVector vs) = SV.unsafeFreeze vs >>= return . BV.convert
+
+-- | concatenate a vector of 'DenseVector's.
+-- If all the dense-vectors are orignally placed consecutively in storage, the result
+-- is simply a type-cast. Otherwise, a new storage is obtained, and dense-vectors are
+-- copied.
+denseVectorConcat :: V.Storable a => BV.Vector (DenseVector a) -> IO (DenseVector a)
+denseVectorConcat vs = do
+  let n = BV.length vs
+      DenseVector (V.MVector sz0 ptr0) = BV.head vs
+  cont <- newIORef True
+  size <- newIORef sz0
+  forM_ [0..n-2] $ \i -> do
+    let DenseVector (V.MVector sz1 ptr1) = vs BV.! i
+        DenseVector (V.MVector sz2 ptr2) = vs BV.! (i+1)
+    modifyIORef cont (&& (V.getPtr ptr1 `advancePtr` sz1) == V.getPtr ptr2)
+    modifyIORef size (+ sz2)
+  cont <- readIORef cont
+  size <- readIORef size
+  if cont
+    then do
+      return $ DenseVector $ V.unsafeFromForeignPtr0 ptr0 size
+    else do
+      nvec@(DenseVector rv) <- newDenseVector size
+      go rv vs
+      return nvec
+  where
+    go vt vs =
+      if BV.null vs
+        then assert (V.length vt == 0) $ return ()
+        else do
+          let DenseVector src = BV.head vs
+              (v1, v2) = V.splitAt (V.length src) vt
+          V.unsafeCopy v1 src
+          go v2 (BV.tail vs)
+
+-- | split a 'DenseVector' into a vector of 'DenseVector's.
+denseVectorSplit :: V.Storable a => Int -> Int -> DenseVector a -> BV.Vector (DenseVector a)
+denseVectorSplit n c (DenseVector v) = assert (V.length v > n * c) $
+  BV.map (\i -> DenseVector (V.unsafeSlice (i*c) c v)) $ BV.enumFromN 0 n
+
+sliceM :: V.Storable a => DenseMatrix a -> (Int, Int) -> DenseVector a
+sliceM (DenseMatrix r c d) (x,y) = assert (x>=0 && x<r && y>=0 && y<c) $ DenseVector v
+  where
+    v = V.unsafeDrop (x*c+y) d
+
+dropV n (DenseVector v) = DenseVector (V.unsafeDrop n v)
+
+copyV (DenseVector v1) (DenseVector v2) len =
+  assert (V.length v1 >= len && V.length v2 >= len) $
+  V.unsafeCopy (V.unsafeTake len v1) (V.unsafeTake len v2)
+
+unsafeReadV :: V.Storable a => DenseVector a -> Int -> IO a
+unsafeReadV (DenseVector v) i = V.unsafeRead v i
+
+unsafeWriteV :: V.Storable a => DenseVector a -> Int -> a -> IO ()
+unsafeWriteV (DenseVector v) i a = V.unsafeWrite v i a
+
+unsafeReadM :: V.Storable a => DenseMatrix a -> (Int, Int) -> IO a
+unsafeReadM (DenseMatrix r c v) (i,j) = assert (i < r && j < c) $ V.unsafeRead v (i*c+j)
+
+unsafeWriteM :: V.Storable a => DenseMatrix a -> (Int, Int) -> a -> IO ()
+unsafeWriteM (DenseMatrix r c v) (i,j) a = assert (i < r && j < c) $ V.unsafeWrite v (i*c+j) a
+
+-- | The Size class provides a interface to tell the dimension of a
+-- dense-vector, dense-matrix, or dense-matrix-array.
+class Size a where
+  type Dim a
+  size :: a -> Dim a
+
+instance V.Storable a => Size (DenseVector a) where
+  type Dim (DenseVector a) = Int
+  size (DenseVector v) = V.length v
+
+instance V.Storable a => Size (DenseMatrix a) where
+  type Dim (DenseMatrix a) = (Int,Int)
+  size (DenseMatrix r c v) = assert (V.length v >= r * c) $ (r,c)
+
+instance V.Storable a => Size (DenseMatrixArray a) where
+  type Dim (DenseMatrixArray a) = (Int,Int,Int)
+  size (DenseMatrixArray n r c v) = assert (V.length v >= n * r * c) $ (n,r,c)
+
+infix 4 :<#, :#>, :<>, :##, :.*, :.+
+infix 0 <<=, <<+
+
+-- | Operations that abstract the low-level details of blas-hs
+data Op :: (* -> *) -> * -> * where
+  -- | vector (as-row) and matrix production
+  (:<#) :: DenseVector a -> DenseMatrix a -> Op DenseVector a
+  -- | matrix and vector (as-column) product
+  (:#>) :: DenseMatrix a -> DenseVector a -> Op DenseVector a
+  -- | matrix and matrix product.
+  -- This is a specially customized matrix matrix product, for the sake of quick
+  -- convolution. The 1st matrix is transposed before multiplication, and the
+  -- result matrix is stored in column-major mode.
+  (:<>) :: DenseMatrix a -> DenseMatrix a -> Op DenseMatrix a
+  -- | vector and vector outer-product
+  (:##) :: DenseVector a -> DenseVector a -> Op DenseMatrix a
+  -- | pairwise product of vector or matrix
+  (:.*) :: c a -> c a -> Op c a
+  -- | pairwise sum of vector or matrix
+  (:.+) :: c a -> c a -> Op c a
+  -- | scale of vector or matrix
+  Scale :: a -> Op c a
+  -- | apply a SIMD-enabled function
+  Apply :: (SIMDPACK a -> SIMDPACK a) -> Op c a
+  -- | zip with a SIMD-enabled function
+  ZipWith :: (SIMDPACK a -> SIMDPACK a -> SIMDPACK a) -> c a -> c a -> Op c a
+  -- | scale the result of some op.
+  -- It is possible to combine scale and many other operations in a single
+  -- BLAS call.
+  Scale' :: a -> Op c a -> Op c a
+  -- | interpret an op to matrix as an op to matrixarray, where each row
+  -- becomes a matrix. This Op is only used internally inside this module
+  UnsafeM2MA :: Op DenseMatrix a -> Op DenseMatrixArray a
+
+-- | Perform an operation
+class AssignTo c a where
+  -- | store the result of a Op to the lhs
+  (<<=) :: c a -> Op c a -> IO ()
+  -- | add the result of a Op to the lhs and store
+  (<<+) :: c a -> Op c a -> IO ()
+
+instance (Numeric a, V.Storable a, SIMDable a) => AssignTo DenseVector a where
+  (DenseVector v) <<= (DenseVector x :<# DenseMatrix r c y) =
+    assert (V.length x == r && V.length v == c) $ gemv_helper Trans r c 1.0 y c x 0.0 v
+
+  (DenseVector v) <<= (DenseMatrix r c x :#> DenseVector y) =
+    assert (V.length y == c && V.length v == r) $ gemv_helper NoTrans r c 1.0 x c y 0.0 v
+
+  (DenseVector v) <<= (DenseVector x :.* DenseVector y) =
+    let sz = V.length v
+    in assert (sz == V.length x && sz == V.length y) $
+       hadamard times v x y
+
+  (DenseVector v) <<= (DenseVector x :.+ DenseVector y) =
+    let sz = V.length v
+    in assert (sz == V.length x && sz == V.length y) $
+       hadamard plus v x y
+
+  (DenseVector v) <<= Scale s =
+    V.unsafeWith v (\pv -> scal (V.length v) s pv 1)
+
+  (DenseVector v) <<= Apply f = foreach f v v
+
+  (DenseVector v) <<= ZipWith f (DenseVector x) (DenseVector y) = hadamard f v x y
+
+  (DenseVector v) <<= Scale' a (DenseMatrix r c x :#> DenseVector y) =
+    assert (V.length y == c && V.length v == r) $ gemv_helper NoTrans r c a x c y 0.0 v
+
+  _ <<= _ = error "Unsupported Op [Vector <<=]."
+
+  (DenseVector v) <<+ (DenseVector x :<# DenseMatrix r c y) =
+    assert (V.length x == r && V.length v == c) $ gemv_helper Trans r c 1.0 y c x 1.0 v
+
+  (DenseVector v) <<+ (DenseMatrix r c x :#> DenseVector y) =
+    assert (V.length y == c && V.length v == r) $ gemv_helper NoTrans r c 1.0 x c y 1.0 v
+
+  (DenseVector v) <<+ Scale' a (DenseMatrix r c x :#> DenseVector y) =
+    assert (V.length y == c && V.length v == r) $ gemv_helper NoTrans r c a x c y 1.0 v
+
+  _ <<+ _ = error "Unsupported Op [Vector <<+]."
+
+instance (Numeric a, V.Storable a, SIMDable a) => AssignTo DenseMatrix a where
+  (DenseMatrix vr vc v) <<= (DenseMatrix xr xc x :<> DenseMatrix yr yc y) =
+    assert (xc == yc && vc == xr && vr == yr) $ do
+      gemm_helper Trans NoTrans xr yr xc 1.0 x xc y xc 0.0 v xr
+
+  (DenseMatrix vr vc v) <<= (DenseMatrix xr xc x :.* DenseMatrix yr yc y) =
+    assert (vr == xr && vr == yr && vc == xc && vc == yc) $ hadamard times v x y
+
+  (DenseMatrix vr vc v) <<= (DenseMatrix xr xc x :.+ DenseMatrix yr yc y) =
+    assert (vr == xr && vr == yr && vc == xc && vc == yc) $ hadamard plus v x y
+
+  (DenseMatrix r c v) <<= Scale s =
+    let sz = V.length v
+    in assert (sz == r * c) $
+       V.unsafeWith v (\pv -> scal sz s pv 1)
+
+  (DenseMatrix r c v) <<= Apply f = (DenseVector v) <<= Apply f
+
+  (DenseMatrix vr vc v) <<= Scale' a (DenseMatrix xr xc x :<> DenseMatrix yr yc y) =
+    assert (xc == yc && vc == xr && vr == yr) $ do
+      gemm_helper Trans NoTrans xr yr xc a x xc y xc 0.0 v xr
+
+  _ <<= _ = error "Unsupported Op [Matrix <<=]."
+
+  (DenseMatrix vr vc v) <<+ (DenseMatrix xr xc x :<> DenseMatrix yr yc y) =
+    assert (xc == yc && vc == xr && vr == yr) $ do
+      gemm_helper Trans NoTrans xr yr xc 1.0 x xc y xc 1.0 v xr
+
+  (DenseMatrix vr vc v) <<+ (DenseVector x :## DenseVector y) =
+    let m = V.length x
+        n = V.length y
+    in assert (m == vr && n == vc) $
+       V.unsafeWith v (\pv ->
+       V.unsafeWith x (\px ->
+       V.unsafeWith y (\py ->
+         geru RowMajor m n 1.0 px 1 py 1 pv n)))
+
+  (DenseMatrix vr vc v)  <<+ Scale' a (DenseMatrix xr xc x :<> DenseMatrix yr yc y) =
+    assert (xc == yc && vc == xr && vr == yr) $ do
+      gemm_helper Trans NoTrans xr yr xc a x xc y xc 1.0 v xr
+
+  _ <<+ _ = error "Unsupported Op [Matrix <<+]."
+
+instance (Numeric a, V.Storable a, SIMDable a) => AssignTo DenseMatrixArray a where
+  ma <<= UnsafeM2MA op = let ma2m (DenseMatrixArray n r c v) = DenseMatrix n (r*c) v
+                         in (ma2m ma) <<= op
+  ma <<= Scale' r (UnsafeM2MA op) = ma <<= UnsafeM2MA (Scale' r op)
+  _ <<= _ = error "Unsupported Op [MatrixArray <<=]."
+  ma <<+ UnsafeM2MA op = let ma2m (DenseMatrixArray n r c v) = DenseMatrix n (r*c) v
+                         in (ma2m ma) <<+ op
+  ma <<+ Scale' r (UnsafeM2MA op) = ma <<+ UnsafeM2MA (Scale' r op)
+  _ <<+ _ = error "Unsupported Op [MatrixArray <<+]."
+
+-- | sum up all elements in the 'DenseMatrix'
+sumElements :: (V.Storable a, Num a) => DenseMatrix a -> IO a
+sumElements (DenseMatrix r c v) = go v (r*c) 0
+  where
+    go v 0  !s = return s
+    go v !n !s = do a <- V.unsafeRead v 0
+                    go (V.unsafeTail v) (n-1) (a+s)
+
+-- | 2D correlation.
+-- Apply a vector of kernels to a dense-matrix with some zero-padding.
+corr2 :: (V.Storable a, Numeric a)
+      => Int                             -- ^ number of 0s padded around
+      -> BV.Vector (DenseMatrix a)       -- ^ vector of kernels
+      -> DenseMatrix a                   -- ^ matrix to be operated
+      -> (Op DenseMatrixArray a -> IO b) -- ^ how to perform the final operation
+      -> IO b
+corr2 p ks m fun = do
+  let k0      = BV.head ks
+      (kr,kc) = size k0
+      (mr,mc) = size m
+      u       = mr - kr + 2*p + 1
+      v       = mc - kc + 2*p + 1
+  zpd <- zero m mr mc p
+  wrk <- newDenseMatrix (u*v) (kr*kc)
+  fill wrk zpd u v kr kc
+  DenseMatrixArray n r c v <- denseMatrixArrayFromVector ks
+  fun $ UnsafeM2MA $ wrk :<> DenseMatrix n (r*c) v
+
+-- | 2D convolution.
+-- Apply a vector of kernels to a dense-matrix with some zero-padding.
+conv2 :: (V.Storable a, Numeric a)
+      => Int                             -- ^ number of 0s padded around
+      -> BV.Vector (DenseMatrix a)       -- ^ vector of kernels
+      -> DenseMatrix a                   -- ^ matrix to be operated
+      -> (Op DenseMatrixArray a -> IO b) -- ^ how to perform the final operation
+      -> IO b
+conv2 p ks m fun = do
+  let k0      = BV.head ks
+      (kr,kc) = size k0
+      (mr,mc) = size m
+      u       = mr - kr + 2*p + 1
+      v       = mc - kc + 2*p + 1
+  zpd <- zero m mr mc p
+  wrk <- newDenseMatrix (u*v) (kr*kc)
+  fill wrk zpd u v kr kc
+  -- copy the kernels, and reverse each.
+  let nk      = BV.length ks
+  knl@(DenseMatrixArray _ _ _ v) <- newDenseMatrixArray nk kr kc
+  forM_ [0..nk-1] $ \i -> do
+    let DenseMatrix _ _ d = denseMatrixArrayAt knl i
+    let DenseMatrix _ _ s = ks BV.! (nk-1-i)
+    V.unsafeCopy d s
+  reverseV v
+  fun $ UnsafeM2MA $ wrk :<> DenseMatrix nk (kr*kc) v
+  where
+    reverseV v = let e = V.length v
+                     m = e `div` 2
+                 in forM_ [0..m] (\i -> V.unsafeSwap v i (e-1-i))
+
+zero m mr mc p = do
+  zpd <- newDenseMatrix (mr+2*p) (mc+2*p)
+  forM_ [0..mr-1] $ \i -> do
+    let t = sliceM zpd (p+i, p)
+        s = sliceM m   (  i, 0)
+    copyV t s mc
+  return zpd
+
+fill wrk@(DenseMatrix _ _ vwrk) m u v kr kc = do
+  refv <- newIORef (DenseVector vwrk)
+  forM_ [0..u-1] $ \i -> do
+    forM_ [0..v-1] $ \j -> do
+      forM_ [0..kr-1] $ \k -> do
+        t <- readIORef refv
+        let s = sliceM m (i+k, j)
+        copyV t s kc
+        writeIORef refv (dropV kc t)
+
+-- | max-pooling, 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 -> DenseMatrix Float -> IO (DenseVector Int, DenseMatrix Float)
+pool 1 mat = do
+  let (r,c) = size mat
+  vi <- newDenseVector (r*c)
+  return (vi, mat)
+pool stride mat = do
+  mxi <- newDenseVector (r'*c')
+  mxv <- newDenseMatrix r' c'
+  forM_ [0..r'-1] $ \i -> do
+    forM_ [0..c'-1] $ \j -> do
+      (n,v) <- unsafeMaxIndEle mat (i*stride) (j*stride) stride stride
+      unsafeWriteV mxi (i*c'+j) n
+      unsafeWriteM mxv (i,j)    v
+  return (mxi,mxv)
+  where
+    (r,c) = size mat
+    r'    = r `div` stride
+    c'    = c `div` stride
+    unsafeMaxIndEle mm x y r c = do
+      mp <- newIORef 0
+      mv <- newIORef (-10000.0)
+      forM_ [0..r-1] $ \ i -> do
+        forM_ [0..c-1] $ \ j -> do
+          v1 <- unsafeReadM mm (x+i, y+j)
+          v0 <- readIORef mv
+          when (v1 > v0) $ do
+            writeIORef mv v1
+            writeIORef mp (i*stride+j)
+      p <- readIORef mp
+      v <- readIORef mv
+      return (p, v)
+
+-- | The reverse of max-pooling.
+unpool :: Int -> DenseVector Int -> DenseMatrix Float -> IO (DenseMatrix Float)
+unpool stride idx mat = do
+  mat' <- newDenseMatrix r' c'
+  forM_ [0..r-1] $ \i -> do
+    forM_ [0..c-1] $ \j -> do
+      pos <- unsafeReadV idx (i*c+j)
+      val <- unsafeReadM mat (i,j)
+      let (oi,oj) = pos `divMod` 2
+      unsafeWriteM mat' (i*stride+oi, j*stride+oj) val
+  return mat'
+  where
+    (r,c) = size mat
+    (r',c') = (r*stride, c*stride)
+
+-- | transpose a vector of 'DenseMatrixArray'
+-- The result is vector of vector of 'DenseMatrix', because the matrices are
+-- no longer placed consecutively in storage.
+transpose :: V.Storable a => BV.Vector (DenseMatrixArray a) -> IO (BV.Vector (BV.Vector (DenseMatrix a)))
+transpose vma = do
+  let DenseMatrixArray n _ _ _  = BV.head vma
+      !vv = BV.map (\i -> BV.map (`denseMatrixArrayAt` i) vma) $ BV.enumFromN 0 n
+  return vv
+
+gemv_helper :: Numeric a
+            => Transpose
+            -> Int -> Int
+            -> a
+            -> V.IOVector a
+            -> Int
+            -> V.IOVector a
+            -> a
+            -> V.IOVector a -> IO ()
+gemv_helper trans row col alpha x lda y beta v =
+  V.unsafeWith x (\px ->
+  V.unsafeWith y (\py ->
+  V.unsafeWith v (\pv ->
+    gemv RowMajor trans row col alpha px lda py 1 beta pv 1)))
+
+gemm_helper :: Numeric a
+            => Transpose
+            -> Transpose
+            -> Int -> Int -> Int
+            -> a
+            -> V.IOVector a
+            -> Int
+            -> V.IOVector a
+            -> Int
+            -> a
+            -> V.IOVector a
+            -> Int
+            -> IO ()
+gemm_helper transA transB rowA colB colA alpha x xlda y ylda beta v vlda =
+  V.unsafeWith x (\px ->
+  V.unsafeWith y (\py ->
+  V.unsafeWith v (\pv -> do
+    gemm ColMajor transA transB rowA colB colA alpha px xlda py ylda beta pv vlda)))
+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,54 @@+import Distribution.Simple
+import Distribution.Simple.Setup
+import Distribution.Simple.LocalBuildInfo
+import Distribution.Simple.Program
+import Distribution.Simple.Program.Run
+import Distribution.Simple.Program.Ar
+import Distribution.Simple.BuildPaths
+import Distribution.System
+import Distribution.Verbosity
+import Distribution.PackageDescription
+import System.FilePath ( (</>) )
+import System.Directory( canonicalizePath, doesDirectoryExist )
+import Control.Monad   ( when, filterM )
+
+main = defaultMainWithHooks simpleUserHooks { buildHook  = myBuildHook }
+
+myBuildHook :: PackageDescription -> LocalBuildInfo -> UserHooks -> BuildFlags -> IO ()
+myBuildHook pkgdesc binfo uh bf = do
+  let buildroot = buildDir binfo
+      vecmod = "compare"
+      sobj = buildroot </> (vecmod ++ ".o")
+      Just dirs = hsSourceDirs . libBuildInfo <$> library pkgdesc
+  dir <- filterM (doesDirectoryExist . (</> "cbits")) dirs
+  if (null dir)
+    then
+      buildHook simpleUserHooks pkgdesc binfo uh bf
+    else do
+      let sdir = head dir
+          ssrc = sdir </> "cbits" </> (vecmod ++ ".ll")
+          verb = fromFlagOrDefault normal (buildVerbosity bf)
+      runProgramInvocation verb $ simpleProgramInvocation "llc" ["-filetype=obj", "-o=" ++ sobj, ssrc]
+      -- generate a ".a" archive for the executable
+      let slib = buildroot </> ("lib" ++ vecmod ++ ".a")
+      createArLibArchive verb binfo slib [sobj]
+      extralib <- canonicalizePath buildroot
+      let pkgdesc' = updatePackageDescription (Nothing, [("t1", libBI)]) pkgdesc
+          libBI    = emptyBuildInfo {extraLibs = [vecmod], extraLibDirs = [extralib]}
+      buildHook simpleUserHooks pkgdesc' binfo uh bf
+      -- however the library is static and doesn't include the vecmod
+      -- we will then explicitly to insert it.
+      let unitId  = componentUnitId (getComponentLocalBuildInfo binfo CLibName)
+          vlibPath = buildroot </> mkLibName     unitId
+          plibPath = buildroot </> mkProfLibName unitId
+          whenVanillaLib = when (withVanillaLib binfo)
+          whenProfLib    = when (withProfLib binfo)
+          Platform hostArch hostOS = hostPlatform binfo
+          args    = case hostOS of
+                       OSX -> ["-q", "-s"]
+                       _   -> ["-q"]
+      (ar, _) <- requireProgram verb arProgram (withPrograms binfo)
+      whenVanillaLib $
+        runProgramInvocation verb $ programInvocation ar (args ++ [vlibPath, sobj])
+      whenProfLib $
+        runProgramInvocation verb $ programInvocation ar (args ++ [plibPath, sobj])
+ Test/Gen.hs view
@@ -0,0 +1,21 @@+{-# LANGUAGE FlexibleInstances, FlexibleContexts #-}
+module Test.Gen where
+import Test.Utils
+import Test.QuickCheck
+import Control.Monad
+import qualified Numeric.LinearAlgebra as L
+import Numeric.LinearAlgebra.Devel
+import qualified Data.Vector.Storable as V
+
+squared_real_matrices :: Int -> Gen (L.Matrix Float)
+squared_real_matrices k = do
+  vs <- sequence (replicate (k*k) arbitrary) :: Gen [Float]
+  return $ L.reshape k $ V.fromList vs
+
+small_matrices :: Gen (L.Matrix Float)
+small_matrices = do
+  n <- choose (2,10)
+  squared_real_matrices n
+
+pair :: Gen a -> Gen b -> Gen (a,b)
+pair = liftM2 (,)
+ Test/S1.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE FlexibleInstances, FlexibleContexts #-}
+module Main where
+import Test.Hspec
+import Test.QuickCheck
+import Numeric.LinearAlgebra
+import qualified Data.Vector.Storable as V
+import Test.Gen
+import Test.Utils
+
+main = hspec $ do
+  describe "Corr Single" $ do
+    it "implements corr2 correct.0" $ do
+      forAll (pair (squared_real_matrices 3) (squared_real_matrices 4)) $
+        \(m1, m2) -> ioProperty $ do
+          r <- test_corr2 2 m1 m2
+          return $ good_corr2 2 m1 m2 `eqShowWhenFail` r
+    it "implements corr2 correct.1" $ do
+      forAll (pair (squared_real_matrices 9) (squared_real_matrices 9)) $
+        \(m1, m2) -> ioProperty $ do
+          r <- test_corr2 4 m1 m2
+          return $ good_corr2 4 m1 m2 `eq` r
+    it "implements corr2 correct.2" $ do
+      forAll (pair (squared_real_matrices 5) (squared_real_matrices 28)) $
+        \(m1, m2) -> ioProperty $ do
+          r <- test_corr2 2 m1 m2
+          return $ eqShowWhenFail (good_corr2 2 m1 m2) r
+    it "implements corr2 correct.3" $ do
+      forAll (pair (choose (2,30)) (pair small_matrices small_matrices)) $
+        \(p, (m1, m2)) -> ioProperty $ do
+          r <- test_corr2 p m1 m2
+          return $ good_corr2 p m1 m2 `eqShowWhenFail` r
+  describe "Corr Many" $ do
+    it "with 2 kernels" $ do
+      forAll (pair (sequence $ replicate 2 $ squared_real_matrices 3) (squared_real_matrices 7)) $
+        \(m1s, m2) -> ioProperty $ do
+          rs <- test_corr2_arr 2 m1s m2
+          return $ conjoin $ zipWith (\m r -> good_corr2 2 m m2 `eqShowWhenFail` r) m1s rs
+    it "with 5 kernels" $ do
+      forAll (pair (sequence $ replicate 4 $ squared_real_matrices 15) (squared_real_matrices 88)) $
+        \(m1s, m2) -> ioProperty $ do
+          rs <- test_corr2_arr 2 m1s m2
+          ss <- return $ map (\m -> good_corr2 2 m m2) m1s
+          return $ conjoin $ zipWith eq rs ss
+
+eqShowWhenFail m1 m2 =
+    whenFail (do let va = flatten m1
+                 let vb = flatten m2
+                 let err x 0 = x
+                     err x y = abs ((x - y) / y)
+                 let ev = (V.zipWith err va vb)
+                     ei = V.maxIndex ev
+                 putStrLn $ "Max error ration: " ++ show (ev V.! ei, va V.! ei, vb V.! ei)
+                 putStrLn $ show m1
+                 putStrLn $ show m2)
+        (m1 `eq` m2)
+ Test/Utils.hs view
@@ -0,0 +1,69 @@+{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
+module Test.Utils where
+import Numeric.LinearAlgebra
+import Control.Exception
+import Control.Monad
+import qualified Data.NeuralNetwork.Backend.BLASHS.Utils as U
+import qualified Numeric.LinearAlgebra as L
+import Numeric.LinearAlgebra.Devel
+import qualified Data.Vector as BV
+import qualified Data.Vector.Storable as V
+import qualified Data.Vector.Storable.Mutable as MV
+import System.IO.Unsafe
+
+asHM (U.DenseMatrix r c v) = L.reshape c $ unsafePerformIO $ V.freeze v
+asDM m = let (r,c) = size m in U.DenseMatrix r c (unsafePerformIO $ V.thaw $ L.flatten m)
+
+good_corr2 :: Int -> L.Matrix Float -> L.Matrix Float -> L.Matrix Float
+
+good_corr2 p k m | w > s     = good_corr2 p m k
+                 | otherwise = corr2 k padded
+  where
+    (w,h) = L.size k
+    (s,t) = L.size m
+    padded = fromBlocks [[z,0,0]
+                        ,[0,m,0]
+                        ,[0,0,z]]
+    z = konst 0 (p, p)
+
+test_corr2 :: Int -> L.Matrix Float -> L.Matrix Float -> IO (L.Matrix Float)
+test_corr2 p k m | w > s     = test_corr2 p m k
+                 | otherwise = do x@(U.DenseMatrixArray _ _ _ vx) <- U.newDenseMatrixArray 1 r c
+                                  k' <- U.DenseMatrix w h <$> V.thaw (flatten k)
+                                  m' <- U.DenseMatrix s t <$> V.thaw (flatten m)
+                                  U.corr2 p (BV.singleton k') m' (x U.<<=)
+                                  reshape c <$> V.freeze vx
+  where
+    (w,h) = L.size k
+    (s,t) = L.size m
+    (r,c) = (s-w+2*p+1, t-h+2*p+1)
+
+test_corr2_arr :: Int -> [L.Matrix Float] -> L.Matrix Float -> IO [L.Matrix Float]
+test_corr2_arr p ks m = do x@(U.DenseMatrixArray _ _ _ vx) <- U.newDenseMatrixArray n r c
+                           print ("test", n, r, c, MV.length vx)
+                           ks' <- mapM (\k -> U.DenseMatrix w h <$> V.thaw (flatten k)) ks
+                           m'  <- U.DenseMatrix s t <$> V.thaw (flatten m)
+                           U.corr2 p (BV.fromList ks') m' (x U.<<=)
+                           let vm = U.denseMatrixArrayToVector x
+                           vhm <- BV.mapM (\(U.DenseMatrix _ _ vx) -> reshape c <$> V.freeze vx) vm
+                           return $ BV.toList vhm
+
+  where
+    n     = length ks
+    (w,h) = L.size (head ks)
+    (s,t) = L.size m
+    (r,c) = (s-w+2*p+1, t-h+2*p+1)
+
+eq :: L.Matrix Float -> L.Matrix Float -> Bool
+eq a b = V.all id $ ratio a b
+
+ratio a b =
+  let va = flatten a
+      vb = flatten b
+      ae :: V.Vector Float
+      ae = V.zipWith (\a b -> abs (a - b)) va vb
+      aa = V.sum ae / fromIntegral (V.length ae)
+      err x 0 = x < 0.1
+      err x y = let e = x-y
+                in (abs (e / y) < 0.02)
+  in V.zipWith err va vb
+ neural-network-blashs.cabal view
@@ -0,0 +1,59 @@+name:                neural-network-blashs
+version:             0.1.0.0
+license-file:        LICENSE
+license:             BSD3
+author:              Jiasen Wu
+maintainer:          jiasenwu@hotmail.com
+homepage:            https://github.com/pierric/neural-network
+bug-reports:         https://github.com/pierric/neural-network/issues
+Category:            AI
+Synopsis:            Yet Another High Performance and Extendable Neural Network in Haskell
+Description:         Provides execution backend of neural network on top of blas-hs.
+Stability:           Experimental
+build-type:          Custom
+cabal-version:       >=1.10
+
+flag vec128
+  Description: Enable 128-bit vector hardware instructions.
+  Default:     False
+  Manual:      True
+flag vec256
+  Description: Enable 256-bit vector hardware instructions.
+  Default:     False
+  Manual:      True
+flag vec512
+  Description: Enable 512-bit vector hardware instructions.
+  Default:     False
+  Manual:      True
+
+library
+  if flag(vec128)
+    hs-source-dirs:     ., vec128
+    ghc-options:        -fllvm
+  else
+    if flag(vec256)
+      hs-source-dirs:   ., vec256
+      ghc-options:      -fllvm -mavx2
+    else
+      if flag(vec512)
+        hs-source-dirs: ., vec512
+        ghc-options:    -fllvm -mavx512
+      else
+        hs-source-dirs: ., novec
+
+  build-depends:       base >= 4.7 && < 5, blas-hs, mwc-random, mtl, vector, constraints, ghc-prim, neural-network-base
+  exposed-modules:     Data.NeuralNetwork.Backend.BLASHS
+  other-modules:       Data.NeuralNetwork.Backend.BLASHS.Layers
+                       Data.NeuralNetwork.Backend.BLASHS.Utils
+                       Data.NeuralNetwork.Backend.BLASHS.SIMD
+  default-language:    Haskell2010
+
+test-suite s1
+  type:               exitcode-stdio-1.0
+  main-is:            Test/S1.hs
+  hs-source-dirs:     ., novec
+  other-modules:      Data.NeuralNetwork.Backend.BLASHS.Utils,
+                      Data.NeuralNetwork.Backend.BLASHS.SIMD
+                      Test.Utils, Test.Gen
+  build-depends:      hspec, QuickCheck, base, hmatrix, vector, blas-hs, neural-network-base
+  default-language:   Haskell2010
+ novec/Data/NeuralNetwork/Backend/BLASHS/SIMD.hs view
@@ -0,0 +1,70 @@+------------------------------------------------------------
+-- |
+-- Module      :  Data.NeuralNetwork.Backend.BLASHS.SIMD
+-- Description :  SIMD based calculations
+-- Copyright   :  (c) 2016 Jiasen Wu
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Jiasen Wu <jiasenwu@hotmail.com>
+-- Stability   :  experimental
+-- Portability :  portable
+--
+--
+-- This module supplies a collection of calculations that
+-- could be implemented on top of SIMD.
+------------------------------------------------------------
+{-# LANGUAGE TypeFamilies, FlexibleContexts #-}
+module Data.NeuralNetwork.Backend.BLASHS.SIMD (
+  SIMDable(..),
+  cost', relu, relu'
+)  where
+
+import Data.Vector.Storable.Mutable as MV
+import Control.Exception
+import qualified Data.NeuralNetwork as B
+
+class SIMDable a where
+  data SIMDPACK a
+  hadamard :: (SIMDPACK a -> SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IOVector a -> IO ()
+  konst    :: a -> SIMDPACK a
+  foreach  :: (SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IO ()
+  plus  :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+  minus :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+  times :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+
+
+instance SIMDable Float where
+  newtype SIMDPACK Float = F { unF :: Float}
+  plus  (F a) (F b) = F (a + b)
+  minus (F a) (F b) = F (a - b)
+  times (F a) (F b) = F (a * b)
+  hadamard op v x y = assert (MV.length x == sz && MV.length y == sz) $ do
+    go sz v x y
+    where
+      sz = MV.length v
+      go 0 _ _ _ = return ()
+      go n z x y = do
+        a <- unsafeRead x 0
+        b <- unsafeRead y 0
+        unsafeWrite z 0 (unF $ op (F a) (F b))
+        go (n-1) (unsafeTail z) (unsafeTail x) (unsafeTail y)
+
+  konst = F
+
+  foreach op v x = assert (sz == MV.length x) $ do
+    go sz v x
+    where
+      sz = MV.length v
+      go 0 _ _ = return ()
+      go n z x = do
+        a <- unsafeRead x 0
+        unsafeWrite z 0 (unF $ op (F a))
+        go (n-1) (unsafeTail z) (unsafeTail x)
+
+-- | SIMD based, RELU and derivative of RELU
+relu, relu' :: SIMDPACK Float -> SIMDPACK Float
+relu  (F a) = F $ B.relu  a
+relu' (F a) = F $ B.relu' a
+
+-- | SIMD based, derivative of error measurement
+cost' :: SIMDPACK Float -> SIMDPACK Float -> SIMDPACK Float
+cost' (F a) (F b) = F (B.cost' a b)
+ vec128/Data/NeuralNetwork/Backend/BLASHS/SIMD.hs view
@@ -0,0 +1,193 @@+------------------------------------------------------------
+-- |
+-- Module      :  Data.NeuralNetwork.Backend.BLASHS.SIMD
+-- Description :  SIMD based calculations
+-- Copyright   :  (c) 2016 Jiasen Wu
+-- License     :  BSD-style (see the file LICENSE)
+-- Maintainer  :  Jiasen Wu <jiasenwu@hotmail.com>
+-- Stability   :  experimental
+-- Portability :  portable
+--
+--
+-- This module supplies a collection of calculations that
+-- could be implemented on top of SIMD.
+------------------------------------------------------------
+{-# LANGUAGE TypeFamilies, FlexibleContexts, FlexibleInstances #-}
+{-# LANGUAGE UnboxedTuples, MagicHash #-}
+{-# LANGUAGE GHCForeignImportPrim, UnliftedFFITypes #-}
+module Data.NeuralNetwork.Backend.BLASHS.SIMD (
+  compareVector,
+  selectVector,
+  SIMDable(..),
+  cost', relu, relu'
+) where
+
+import Data.Vector.Storable.Mutable as MV
+import qualified Data.Vector.Storable as SV
+import Control.Exception
+import Control.Monad
+
+import GHC.Prim
+import GHC.Base
+import GHC.Exts
+import GHC.Ptr (Ptr(..))
+import Foreign.Storable (Storable(..))
+
+foreign import prim "vfcomp_oge" fcomp_oge :: FloatX4# -> FloatX4# -> Word32X4#
+foreign import prim "vselect"    select    :: Word32X4# -> FloatX4# -> FloatX4# -> FloatX4#
+
+data Word32X4 = Word32X4 Word32X4#
+
+data CompareFunc = GE
+
+class SIMDVector v => Comparable v where
+  compareVector :: CompareFunc -> v -> v -> Word32X4
+  selectVector  :: Word32X4 -> v -> v -> v
+
+instance Comparable (SIMDPACK Float) where
+  compareVector GE (FloatX4 x) (FloatX4 y) = Word32X4 (fcomp_oge x y)
+  selectVector (Word32X4 s) (FloatX4 x) (FloatX4 y) = FloatX4 (select s x y)
+
+class SIMDable a where
+  data SIMDPACK a
+  hadamard :: (SIMDPACK a -> SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IOVector a -> IO ()
+  konst    :: a -> SIMDPACK a
+  foreach  :: (SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IO ()
+  plus     :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+  minus    :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+  times    :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+
+instance SIMDable Float where
+  data SIMDPACK Float = FloatX4 FloatX4#
+  plus   (FloatX4 a) (FloatX4 b) = FloatX4 (plusFloatX4#  a b)
+  minus  (FloatX4 a) (FloatX4 b) = FloatX4 (minusFloatX4# a b)
+  times  (FloatX4 a) (FloatX4 b) = FloatX4 (timesFloatX4# a b)
+  hadamard op v x y = assert (MV.length x == sz && MV.length y == sz) $ do
+    let sv = unsafeCast v :: IOVector (SIMDPACK Float)
+        sx = unsafeCast x :: IOVector (SIMDPACK Float)
+        sy = unsafeCast y :: IOVector (SIMDPACK Float)
+    go (MV.length sv) sv sx sy
+    let rm = sz `mod` 4
+        rn = sz - rm
+        rv = unsafeDrop rn v
+        rx = unsafeDrop rn x
+        ry = unsafeDrop rn y
+    when (rm /= 0) $ rest rm rv rx ry
+    where
+      sz = MV.length v
+      go 0 _ _ _ = return ()
+      go n z x y = do
+        a <- unsafeRead x 0
+        b <- unsafeRead y 0
+        unsafeWrite z 0 (op a b)
+        go (n-1) (unsafeTail z) (unsafeTail x) (unsafeTail y)
+      rest n z x y = do
+        sx <- SV.unsafeFreeze x
+        sy <- SV.unsafeFreeze y
+        let vx = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sx
+            vy = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sy
+            (vz0,vz1,vz2,_) = unpackVector (op vx vy)
+        unsafeWrite z 0 vz0
+        when (n > 1) $ do
+          unsafeWrite z 1 vz1
+          when (n > 2) $ do
+            unsafeWrite z 2 vz2
+
+  konst = broadcastVector
+
+  foreach op v x = assert (sz == MV.length x) $ do
+    let sv = unsafeCast v :: IOVector (SIMDPACK Float)
+        sx = unsafeCast x :: IOVector (SIMDPACK Float)
+    go (MV.length sv) sv sx
+    let rm = sz `mod` 4
+        rn = sz - rm
+        rv = unsafeDrop rn v
+        rx = unsafeDrop rn x
+    when (rm /= 0) $ rest rm rv rx
+    where
+      sz = MV.length v
+      go 0 _ _ = return ()
+      go n z x = do
+        a <- unsafeRead x 0
+        unsafeWrite z 0 (op a)
+        go (n-1) (unsafeTail z) (unsafeTail x)
+      rest n z x = do
+        sx <- SV.unsafeFreeze x
+        let vx = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sx
+            (vz0,vz1,vz2,_) = unpackVector (op vx)
+        unsafeWrite z 0 vz0
+        when (n > 1) $ do
+          unsafeWrite z 1 vz1
+          when (n > 2) $ do
+            unsafeWrite z 2 vz2
+
+-- | SIMD based, RELU and derivative of RELU
+relu, relu' :: SIMDPACK Float -> SIMDPACK Float
+relu  x = let v0 = broadcastVector 0
+          in selectVector (compareVector GE x v0) x v0
+relu' x = let v0 = broadcastVector 0
+              v1 = broadcastVector 1
+          in selectVector (compareVector GE v0 x) v0 v1
+
+-- | SIMD based, derivative of error measurement
+cost' :: SIMDPACK Float -> SIMDPACK Float -> SIMDPACK Float
+cost' a y = selectVector (compareVector GE a y)
+              (selectVector (compareVector GE y (broadcastVector 1))
+                (broadcastVector 0)
+                (minus a y))
+              (minus a y)
+
+instance Storable (SIMDPACK Float) where
+    sizeOf x     = vectorSize x * elementSize x
+    alignment    = sizeOf
+    peek (Ptr a) = IO $ \s -> let (# s', r #) = readFloatX4OffAddr# a 0# s in (# s', FloatX4 r #)
+    poke (Ptr a) (FloatX4 b) = IO $ \s -> (# writeFloatX4OffAddr# a 0# b s, () #)
+
+class SIMDVector v where
+    -- | Type of the elements in the vector
+    type Elem v
+    -- | Type used to pack or unpack the vector
+    type ElemTuple v
+    -- | Vector with all elements initialized to zero.
+    nullVector       :: v
+    -- | Number of components (scalar elements) in the vector. The argument is not evaluated.
+    vectorSize       :: v -> Int
+    -- | Size of each (scalar) element in the vector in bytes. The argument is not evaluated.
+    elementSize      :: v -> Int
+    -- | Broadcast a scalar to all elements of a vector.
+    broadcastVector  :: Elem v -> v
+    -- | Insert a scalar at the given position (starting from 0) in a vector. If the index is outside of the range an exception is thrown.
+    insertVector     :: v -> Elem v -> Int -> v
+    insertVector v e i | i < 0            = error $ "insertVector: negative argument: " ++ show i
+                       | i < vectorSize v = unsafeInsertVector v e i
+                       | otherwise        = error $ "insertVector: argument too large: " ++ show i
+    -- | Insert a scalar at the given position (starting from 0) in a vector. If the index is outside of the range the behavior is undefined.
+    unsafeInsertVector     :: v -> Elem v -> Int -> v
+    -- | Pack some elements to a vector.
+    packVector       :: ElemTuple v -> v
+    -- | Unpack a vector.
+    unpackVector     :: v -> ElemTuple v
+
+instance SIMDVector (SIMDPACK Float) where
+    type Elem (SIMDPACK Float) = Float
+    type ElemTuple (SIMDPACK Float) = (Float, Float, Float, Float)
+    nullVector         = broadcastVector 0
+    vectorSize  _      = 4
+    elementSize _      = 4
+    broadcastVector    = broadcastFloatX4
+    unsafeInsertVector = unsafeInsertFloatX4
+    packVector         = packFloatX4
+    unpackVector       = unpackFloatX4
+
+{-# INLINE broadcastFloatX4 #-}
+broadcastFloatX4 (F# x) = FloatX4 (broadcastFloatX4# x)
+
+{-# INLINE packFloatX4 #-}
+packFloatX4 (F# x1, F# x2, F# x3, F# x4) = FloatX4 (packFloatX4# (# x1, x2, x3, x4 #))
+
+{-# INLINE unpackFloatX4 #-}
+unpackFloatX4 (FloatX4 m1) = case unpackFloatX4# m1 of
+    (# x1, x2, x3, x4 #) -> (F# x1, F# x2, F# x3, F# x4)
+
+{-# INLINE unsafeInsertFloatX4 #-}
+unsafeInsertFloatX4 (FloatX4 m1) (F# y) _i@(I# ip) = FloatX4 (insertFloatX4# m1 y (ip -# 0#))
+ vec256/Data/NeuralNetwork/Backend/BLASHS/SIMD.hs view
@@ -0,0 +1,164 @@+{-# LANGUAGE TypeFamilies, FlexibleContexts, FlexibleInstances #-}
+{-# LANGUAGE UnboxedTuples, MagicHash #-}
+{-# LANGUAGE GHCForeignImportPrim, UnliftedFFITypes #-}
+module Data.NeuralNetwork.Backend.BLASHS.SIMD where
+
+import Data.Vector.Storable.Mutable as MV
+import qualified Data.Vector.Storable as SV
+import Control.Exception
+import Control.Monad
+
+import GHC.Prim
+import GHC.Base
+import GHC.Exts
+import GHC.Ptr (Ptr(..))
+import Foreign.Storable (Storable(..))
+
+foreign import prim "vfcomp_oge" fcomp_oge :: FloatX8# -> FloatX8# -> Word32X8#
+foreign import prim "vselect"    select    :: Word32X8# -> FloatX8# -> FloatX8# -> FloatX8#
+
+data Word32X8 = Word32X8 Word32X8#
+
+data CompareFunc = GE
+
+class SIMDVector v => Comparable v where
+  compareVector :: CompareFunc -> v -> v -> Word32X8
+  selectVector  :: Word32X8 -> v -> v -> v
+
+instance Comparable (SIMDPACK Float) where
+  compareVector GE (FloatX8 x) (FloatX8 y) = Word32X8 (fcomp_oge x y)
+  selectVector (Word32X8 s) (FloatX8 x) (FloatX8 y) = FloatX8 (select s x y)
+
+class SIMDable a where
+  data SIMDPACK a
+  hadamard :: (SIMDPACK a -> SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IOVector a -> IO ()
+  konst    :: a -> SIMDPACK a
+  foreach  :: (SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IO ()
+  plus     :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+  minus    :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+  times    :: SIMDPACK a -> SIMDPACK a -> SIMDPACK a
+
+instance SIMDable Float where
+  data SIMDPACK Float = FloatX8 FloatX8#
+  plus   (FloatX8 a) (FloatX8 b) = FloatX8 (plusFloatX8#  a b)
+  minus  (FloatX8 a) (FloatX8 b) = FloatX8 (minusFloatX8# a b)
+  times  (FloatX8 a) (FloatX8 b) = FloatX8 (timesFloatX8# a b)
+  hadamard op v x y = assert (MV.length x == sz && MV.length y == sz) $ do
+    let sv = unsafeCast v :: IOVector (SIMDPACK Float)
+        sx = unsafeCast x :: IOVector (SIMDPACK Float)
+        sy = unsafeCast y :: IOVector (SIMDPACK Float)
+    go (MV.length sv) sv sx sy
+    let rm = sz `mod` 8
+        rn = sz - rm
+        rv = unsafeDrop rn v
+        rx = unsafeDrop rn x
+        ry = unsafeDrop rn y
+    when (rm /= 0) $ rest rm rv rx ry
+    where
+      sz = MV.length v
+      go 0 _ _ _ = return ()
+      go n z x y = do
+        a <- unsafeRead x 0
+        b <- unsafeRead y 0
+        unsafeWrite z 0 (op a b)
+        go (n-1) (unsafeTail z) (unsafeTail x) (unsafeTail y)
+      rest n z x y = do
+        sx <- SV.unsafeFreeze x
+        sy <- SV.unsafeFreeze y
+        let vx = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sx
+            vy = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sy
+            (vz0,vz1,vz2,vz3,vz4,vz5,vz6,_) = unpackVector (op vx vy)
+        forM_ (zip [0..n-1] [vz0,vz1,vz2,vz3,vz4,vz5,vz6]) $ uncurry (unsafeWrite z)
+
+  konst = broadcastVector
+
+  foreach op v x = assert (sz == MV.length x) $ do
+    let sv = unsafeCast v :: IOVector (SIMDPACK Float)
+        sx = unsafeCast x :: IOVector (SIMDPACK Float)
+    go (MV.length sv) sv sx
+    let rm = sz `mod` 8
+        rn = sz - rm
+        rv = unsafeDrop rn v
+        rx = unsafeDrop rn x
+    when (rm /= 0) $ rest rm rv rx
+    where
+      sz = MV.length v
+      go 0 _ _ = return ()
+      go n z x = do
+        a <- unsafeRead x 0
+        unsafeWrite z 0 (op a)
+        go (n-1) (unsafeTail z) (unsafeTail x)
+      rest n z x = do
+        sx <- SV.unsafeFreeze x
+        let vx = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sx
+            (vz0,vz1,vz2,vz3,vz4,vz5,vz6,_) = unpackVector (op vx)
+        forM_ (zip [0..n-1] [vz0,vz1,vz2,vz3,vz4,vz5,vz6]) $ uncurry (unsafeWrite z)
+
+relu, relu' :: SIMDPACK Float -> SIMDPACK Float
+relu  x = let v0 = broadcastVector 0
+          in selectVector (compareVector GE x v0) x v0
+relu' x = let v0 = broadcastVector 0
+              v1 = broadcastVector 1
+          in selectVector (compareVector GE v0 x) v0 v1
+
+cost' :: SIMDPACK Float -> SIMDPACK Float -> SIMDPACK Float
+cost' a y = selectVector (compareVector GE a y)
+              (selectVector (compareVector GE y (broadcastVector 1))
+                (broadcastVector 0)
+                (minus a y))
+              (minus a y)
+
+instance Storable (SIMDPACK Float) where
+    sizeOf x     = vectorSize x * elementSize x
+    alignment    = sizeOf
+    peek (Ptr a) = IO $ \s -> let (# s', r #) = readFloatX8OffAddr# a 0# s in (# s', FloatX8 r #)
+    poke (Ptr a) (FloatX8 b) = IO $ \s -> (# writeFloatX8OffAddr# a 0# b s, () #)
+
+class SIMDVector v where
+    -- | Type of the elements in the vector
+    type Elem v
+    -- | Type used to pack or unpack the vector
+    type ElemTuple v
+    -- | Vector with all elements initialized to zero.
+    nullVector       :: v
+    -- | Number of components (scalar elements) in the vector. The argument is not evaluated.
+    vectorSize       :: v -> Int
+    -- | Size of each (scalar) element in the vector in bytes. The argument is not evaluated.
+    elementSize      :: v -> Int
+    -- | Broadcast a scalar to all elements of a vector.
+    broadcastVector  :: Elem v -> v
+    -- | Insert a scalar at the given position (starting from 0) in a vector. If the index is outside of the range an exception is thrown.
+    insertVector     :: v -> Elem v -> Int -> v
+    insertVector v e i | i < 0            = error $ "insertVector: negative argument: " ++ show i
+                       | i < vectorSize v = unsafeInsertVector v e i
+                       | otherwise        = error $ "insertVector: argument too large: " ++ show i
+    -- | Insert a scalar at the given position (starting from 0) in a vector. If the index is outside of the range the behavior is undefined.
+    unsafeInsertVector     :: v -> Elem v -> Int -> v
+    -- | Pack some elements to a vector.
+    packVector       :: ElemTuple v -> v
+    -- | Unpack a vector.
+    unpackVector     :: v -> ElemTuple v
+
+instance SIMDVector (SIMDPACK Float) where
+    type Elem (SIMDPACK Float) = Float
+    type ElemTuple (SIMDPACK Float) = (Float, Float, Float, Float, Float, Float, Float, Float)
+    nullVector         = broadcastVector 0
+    vectorSize  _      = 8
+    elementSize _      = 8
+    broadcastVector    = broadcastFloatX8
+    unsafeInsertVector = unsafeInsertFloatX8
+    packVector         = packFloatX8
+    unpackVector       = unpackFloatX8
+
+{-# INLINE broadcastFloatX8 #-}
+broadcastFloatX8 (F# x) = FloatX8 (broadcastFloatX8# x)
+
+{-# INLINE packFloatX8 #-}
+packFloatX8 (F# x1, F# x2, F# x3, F# x4, F# x5, F# x6, F# x7, F# x8 ) = FloatX8 (packFloatX8# (# x1, x2, x3, x4, x5, x6, x7, x8 #))
+
+{-# INLINE unpackFloatX8 #-}
+unpackFloatX8 (FloatX8 m1) = case unpackFloatX8# m1 of
+    (# x1, x2, x3, x4, x5, x6, x7, x8 #) -> (F# x1, F# x2, F# x3, F# x4, F# x5, F# x6, F# x7, F# x8)
+
+{-# INLINE unsafeInsertFloatX8 #-}
+unsafeInsertFloatX8 (FloatX8 m1) (F# y) _i@(I# ip) = FloatX8 (insertFloatX8# m1 y (ip -# 0#))
+ vec512/Data/NeuralNetwork/Backend/BLASHS/SIMD.hs view
@@ -0,0 +1,66 @@+{-# LANGUAGE TypeFamilies, FlexibleContexts, BangPatterns #-}
+module Data.NeuralNetwork.Backend.BLASHS.SIMD where
+
+import Data.Vector.Storable.Mutable as MV
+import qualified Data.Vector.Storable as SV
+import Data.Primitive.SIMD
+import Control.Exception
+import Control.Monad
+
+class Num (SIMDPACK a) => SIMDable a where
+  type SIMDPACK a
+  konst :: a -> SIMDPACK a
+  foreach  :: (SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IO ()
+  hadamard :: (SIMDPACK a -> SIMDPACK a -> SIMDPACK a) -> IOVector a -> IOVector a -> IOVector a -> IO ()
+
+instance SIMDable Float where
+  type SIMDPACK Float = FloatX16
+  hadamard op v x y = assert (MV.length x == sz && MV.length y == sz) $ do
+    let sv = unsafeCast v :: IOVector FloatX16
+        sx = unsafeCast x :: IOVector FloatX16
+        sy = unsafeCast y :: IOVector FloatX16
+    go (MV.length sv) sv sx sy
+    let rm = sz `mod` 16
+        rn = sz - rm
+        rv = unsafeDrop rn v
+        rx = unsafeDrop rn x
+        ry = unsafeDrop rn y
+    when (rm /= 0) $ rest rm rv rx ry
+    where
+      sz = MV.length v
+      go 0 _ _ _ = return ()
+      go !n !z !x !y = do
+        a <- unsafeRead x 0
+        b <- unsafeRead y 0
+        unsafeWrite z 0 (op a b)
+        go (n-1) (unsafeTail z) (unsafeTail x) (unsafeTail y)
+      rest n z x y = do
+        sx <- SV.unsafeFreeze x
+        sy <- SV.unsafeFreeze y
+        let vx = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sx
+            vy = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sy
+            (vz0,vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,vz9,vzA,vzB,vzC,vzD,vzE,_) = unpackVector (op vx vy)
+        forM_ (zip [0..n-1] [vz0,vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,vz9,vzA,vzB,vzC,vzD,vzE]) $ uncurry (unsafeWrite z)
+
+  foreach op v x = assert (MV.length v == MV.length x) $ do
+    let sv = unsafeCast v :: IOVector FloatX16
+        sx = unsafeCast x :: IOVector FloatX16
+    go (MV.length sv) sv sx
+    let rm = sz `mod` 16
+        rn = sz - rm
+        rv = unsafeDrop rn v
+        rx = unsafeDrop rn x
+    when (rm /= 0) $ rest rm rv rx
+    where
+      sz = MV.length v
+      go 0 _ _ = return ()
+      go !n !z !x = do
+        a <- unsafeRead x 0
+        unsafeWrite z 0 (op a)
+        go (n-1) (unsafeTail z) (unsafeTail x)
+      rest n z x = do
+        sx <- SV.unsafeFreeze x
+        let vx = SV.ifoldl' (\v i a -> unsafeInsertVector v a i) nullVector sx
+            (vz0,vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,vz9,vzA,vzB,vzC,vzD,vzE,_) = unpackVector (op vx)
+        forM_ (zip [0..n-1] [vz0,vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,vz9,vzA,vzB,vzC,vzD,vzE]) $ uncurry (unsafeWrite z)
+  konst = broadcastVector