diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Marco Zocca (c) 2017
+
+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 Marco Zocca nor the names of other
+      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
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,12 @@
+# taco
+
+Haskell implementation of the Tensor Algebra COmpiler (`taco`) [1]
+
+**NB: Experimental and not yet ready for general use**
+
+
+
+
+## References
+
+[1] F. Kjolstad et al., Proc. ACM Program. Lang., Vol.1, OOPSLA, Article 77, October 2017
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/src/Data/Dim.hs b/src/Data/Dim.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Dim.hs
@@ -0,0 +1,28 @@
+{-# language GADTs, TypeOperators #-}
+module Data.Dim where
+
+import Data.Vector.Unboxed as V
+
+
+-- * Dimension metadata
+
+-- | To define a /dense/ dimension we only need the dimensionality parameter
+newtype Dd i = Dd { dDim :: i } deriving (Eq, Show)
+
+
+-- | To define a /sparse/ dimension we need a cumulative array, an index array and a dimensionality parameter
+data Sd i = Sd {
+      sCml :: Maybe (V.Vector i)
+    , sIdx :: V.Vector i
+    , sDim :: i }
+  deriving (Eq, Show)
+
+-- | A tensor dimension can be either dense or sparse.
+--
+-- Example: the CSR format is /dense/ in the first index (rows) and /sparse/ in the second index (columns)
+
+-- newtype DMD i = DMD (Either (DMDDense i) (DMDSparse i)) deriving (Eq, Show)
+
+
+  
+  
diff --git a/src/Data/Shape.hs b/src/Data/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Shape.hs
@@ -0,0 +1,102 @@
+{-# language GADTs, TypeOperators #-}
+
+module Data.Shape where
+
+import Data.Monoid
+import Data.Int (Int32)
+-- import GHC.TypeLits
+-- import GHC.Natural
+
+import qualified Data.Vector.Unboxed as VU
+
+import qualified Data.Dim as Dim
+
+
+data Z
+data sh :# e -- dense
+data sh :. e -- sparse
+
+-- | A statically-typed tensor shape parameter that supports both sparse and dense dimensions.
+-- Dimensions are indexed with Int32 indices, which should be enough for most applications.
+data Sh sh where 
+  Z :: Sh Z
+  -- | Constructor for a dense dimension  
+  D :: Sh sh -> Dim.Dd Int32 -> Sh (sh :# Int32)
+  -- | Constructor for a sparse dimension 
+  S :: Sh sh -> Dim.Sd Int32 -> Sh (sh :. Int32) 
+
+type D1 = Z :# Int32
+type D2 = (Z :# Int32) :# Int32
+type CSR = (Z :# Int32) :. Int32
+type COO = (Z :. Int32) :. Int32
+
+
+instance Show (Sh sh) where
+  show Z = ""
+  show (D sh (Dim.Dd m)) = unwords [show m, show sh]
+  show (S sh (Dim.Sd _ ix n)) = showSparse ix n <> show sh where
+    showSparse ixx nn = show (VU.length ixx, nn)
+
+instance Eq (Sh sh) where
+  Z == Z = True
+  (sh `D` d) == (sh2 `D` d2) = d == d2 && (sh == sh2)
+  (sh `S` s) == (sh2 `S` s2) = s == s2 && (sh == sh2)
+
+-- | Rank of a shape (i.e. number of dimensions)
+rank :: Sh sh -> Int
+rank Z = 0
+rank (D sh _) = 1 + rank sh
+rank (S sh _) = 1 + rank sh
+
+-- | Dimension of a shape (i.e. list of dimension sizes)
+dim :: Sh sh -> [Integer]
+dim Z = []
+dim (D sh (Dim.Dd m)) = toInteger m : dim sh
+dim (S sh (Dim.Sd _ _ m)) = toInteger m : dim sh
+
+
+
+
+  
+-- | Shape of a dense rank-2 tensor (a matrix)
+mkD2 :: Int32 -> Int32 -> Sh D2
+mkD2 m n = (Z `D` Dim.Dd m) `D` Dim.Dd n
+
+-- | Shape of a rank-2 CSR tensor (dense in the first index, sparse in the second)
+mkCSR :: Int32 -> Int32 -> VU.Vector Int32 -> VU.Vector Int32 -> Sh CSR
+mkCSR m n icml iidx = (Z `D` Dim.Dd m) `S` Dim.Sd (Just icml) iidx n
+
+-- | Shape of a rank-2 COO tensor (sparse in both indices)
+mkCOO :: Int32 -> Int32 -> VU.Vector Int32 -> VU.Vector Int32 -> Sh COO
+mkCOO m n vi vj = (Z `S` Dim.Sd Nothing vi m) `S` Dim.Sd Nothing vj n
+
+
+
+
+
+-- | Playground
+
+-- -- | An index of dimension zero
+-- data Z  = Z
+--         deriving (Show, Read, Eq, Ord)
+
+-- -- | Our index type, used for both shapes and indices.
+-- infixl 3 :.
+-- data tail :. head = !tail :. !head
+--         deriving (Show, Read, Eq, Ord)
+
+
+
+-- class Eq sh => Shape sh where
+--   rank :: sh -> Int
+--   -- listOfShape :: sh -> [Int]
+--   -- shapeOfList :: [Int] -> sh
+
+-- instance (Eq i, Shape sh) => Shape (sh :. i) where
+--   rank (sh :. _) = rank sh + 1
+
+-- instance Shape Z where
+--   rank _ = 0
+--   -- listOfShape _ = []
+--   -- shapeOfList [] = Z
+--   -- shapeOfList _ = error $ stage ++ ".fromList: non-empty list when converting to Z."
diff --git a/src/Data/Tensor.hs b/src/Data/Tensor.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor.hs
@@ -0,0 +1,152 @@
+{-# language GADTs #-}
+{-# language DeriveFunctor #-}
+{-# language TypeOperators #-}
+module Data.Tensor (
+  -- * Tensor type
+  Tensor(..), shape, nnz,
+  -- * Shape type
+  Sh(..),
+  -- * Dimension types
+  Dim.Dd(..), Dim.Sd(..)) where
+
+import qualified Data.Vector as V
+-- import qualified Data.Vector.Unboxed as VU
+
+-- import Data.Word (Word32, Word64)
+-- import Data.Int (Int32)
+
+import Data.Shape (Sh(..), dim, rank,
+                   Z,
+                   D1, D2, CSR, COO, mkD2, mkCSR, mkCOO)
+import qualified Data.Dim as Dim
+
+
+
+{- |
+IN: Tensor reduction syntax (Einstein notation)
+
+OUT: stride program (how to read/write memory)
+
+
+taco compiles a tensor expression (e.g. C = A_{ijk}B_{k} ) into a series of nested loops.
+
+dimensions : can be either dense or sparse
+
+internally, tensor data is stored in /dense/ vectors
+
+"contract A_{ijk}B_{k} over the third index"
+
+-}
+
+
+
+-- | A generic tensor type, polymorphic in the container type as well
+data GTensor c i a where
+  GTensor :: Sh i -> c a -> GTensor c (Sh i) a
+  
+mkGT :: Sh i -> c a -> GTensor c (Sh i) a
+mkGT = GTensor
+
+
+-- | The 'Tensor' type. Tensor data entries are stored as one single array
+data Tensor i a where
+  T :: Sh i -> V.Vector a -> Tensor (Sh i) a 
+
+-- | Construct a tensor given a shape and a vector of entries
+mkT :: Sh i -> V.Vector a -> Tensor (Sh i) a
+mkT = T
+
+instance Functor (Tensor i) where
+  fmap f (T sh v) = T sh (f <$> v)
+
+-- liftA2' :: (a -> a -> b) -> Tensor i a -> Tensor i a -> Tensor i a 
+-- liftA2' f (T sh1 v1) (T sh2 v2) = mkT sh1 (V.zipWith f v1 v2)
+
+pure' :: a -> Tensor (Sh Z) a
+pure' = mkT Z . V.singleton
+
+instance (Eq a) => Eq (Tensor i a) where
+  (T sh1 d1) == (T sh2 d2) = sh1 == sh2 && d1 == d2
+
+instance (Show a) => Show (Tensor i a) where
+  show (T sh d) = unwords [show sh, show $ V.take 5 d, "..."]
+
+
+ 
+
+-- | Access the shape of a tensor
+shape :: Tensor sh a -> sh
+shape (T sh _) = sh
+
+-- | Number of nonzero tensor elements
+nnz :: Tensor i a -> Int
+nnz (T _ td) = V.length td
+
+
+
+
+
+
+-- * A possible abstract syntax
+
+-- data Index i where
+--   I1 :: i -> Index i
+--   I2 :: i -> i -> Index (i, i)
+
+-- -- | Expressions with tensor operands, e.g. "contract A_{ijk}B_{k} over the third index"
+
+-- -- mkConstE = Const <$> mkT
+
+-- data Expr a where
+--   Const :: Tensor (Sh i) a -> Expr (Tensor (Sh i) a)
+
+-- data Expr i a where
+--   -- Const :: a -> Expr a
+--   Contract :: Index i -> Expr (Sh i) a -> Expr (Sh i) a -> Expr (Sh i) a
+--   -- (:*:) :: Expr a -> Expr a -> Expr a
+--   -- (:+:) :: Expr a -> Expr a -> Expr a
+
+-- eval (Const x) = x
+-- eval (Contract ixs a b) = undefined
+
+
+
+-- data Expr a =
+--     Const a
+--   | Contract Int (Expr a) (Expr a)
+--   -- | Expr a :+: Expr a
+--   -- | Expr a :*: Expr a
+--   -- | Expr a :-: Expr a
+--   -- | Expr a :/: Expr a
+--   deriving (Eq, Show)
+
+-- -- | trivial recursive evaluation function
+-- eval :: Num t => Expr t -> t
+-- eval (Const x) = x
+-- eval (a :+: b) = eval a + eval b
+-- eval (a :*: b) = eval a * eval b
+
+
+
+
+-- | GADT syntax
+
+-- data Expr a where
+--   Const :: a -> Expr a 
+--   -- ^ Sum (elementwise) two expressions
+--   (:+:) :: Expr a -> Expr a -> Expr a
+--   -- ^ Multiply (elementwise) two expressions
+--   (:*:) :: Expr a -> Expr a -> Expr a
+--   -- ^ Subtract (elementwise) two expressions
+--   (:-:) :: Expr a -> Expr a -> Expr a
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/Data/TensorTest.hs b/src/Data/TensorTest.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/TensorTest.hs
@@ -0,0 +1,148 @@
+{-# language MultiParamTypeClasses, TypeFamilies, FlexibleContexts #-}
+{-# language TypeOperators #-}
+{-# language GADTs #-}
+module Data.TensorTest where
+
+import Data.List (splitAt, unfoldr)
+
+import qualified Data.Vector.Unboxed as V
+-- import Data.Tensor
+
+
+data Exp a where
+  Lift :: (a -> b) -> Exp a -> Exp b
+
+
+-- | tensor elements may be indexed
+class Ord (Ix a) => Elem a where
+  type Ix a :: *
+  type Ev a :: *
+  proj :: a -> (Ix a, Ev a)
+  prod :: Ix a -> Ev a -> a
+
+
+data E1 a = E1 Int a deriving Show
+
+instance Elem (E1 a) where
+  type Ix (E1 a) = Int
+  type Ev (E1 a) = a
+  proj (E1 i x) = (i, x)
+  prod i x = E1 i x
+
+data E2 a = E2 Int Int a deriving Show
+
+instance Elem (E2 a) where
+  type Ix (E2 a) = (Int, Int)
+  type Ev (E2 a) = a
+  proj (E2 i j x) = ((i, j), x)
+  prod (i, j) x = E2 i j x
+
+
+
+
+      
+
+
+
+
+spUnion' :: Elem a => (Ev a -> Ev a -> Ev a) -> [a] -> [a] -> [a]
+spUnion' ff = go where
+  go [] y = y
+  go x [] = x
+  go xv@(elx:xs) yv@(ely:ys) =
+    let
+      (ix, elvx) = proj elx
+      (iy, elvy) = proj ely
+    in 
+      case compare ix iy of
+        EQ -> (prod ix (ff elvx elvy)) : go xs ys
+        LT -> prod ix elvx             : go xs yv
+        GT -> prod iy elvy             : go xv ys
+        
+
+
+
+
+-- | sparse vector /union/ as 2-way merge
+spUnion :: Ord i => (a -> a -> a) -> [(i, a)] -> [(i, a)] -> [(i, a)]
+spUnion ff = go where
+  go [] y = y
+  go x [] = x
+  go xv@((i,x):xs) yv@((j,y):ys) = 
+    case compare i j of EQ -> (i, ff x y) : go xs ys
+                        LT -> (i, x)      : go xs yv
+                        GT -> (j, y)      : go xv ys
+
+-- | sparse vector /component-wise intersection/ as 2-way merge
+spIntersect :: Ord i => (a -> a -> a) -> [(i, a)] -> [(i, a)] -> [(i, a)]
+spIntersect gg = go where
+  go [] _  = []
+  go _  [] = []
+  go xv@((i,x):xs) yv@((j,y):ys) =
+    case compare i j of EQ -> (i, gg x y) : go xs ys
+                        LT -> go xs yv
+                        GT -> go xv ys
+
+spAdd :: (Num a, Ord i) => [(i, a)] -> [(i, a)] -> [(i, a)]  
+spAdd = spUnion (+)
+
+spMul :: (Num a, Ord i) => [(i, a)] -> [(i, a)] -> [(i, a)]
+spMul = spIntersect (*)
+
+
+v0, v1 :: [(Int, Int)]
+v0 = [(0, 1), (2, 2), (5, 1)]
+v1 = [(0, 2), (1, 3), (2, 3), (4, 1)]
+
+
+
+-- | "AND"
+conjunction :: (a -> b -> c) -> Maybe a -> Maybe b -> Maybe c
+conjunction f a b = case (a, b) of
+  (Just x, Just y)  -> Just $ f x y
+  _                 -> Nothing
+
+-- | "OR"
+disjunction :: (t -> t -> t) -> Maybe t -> Maybe t -> Maybe t
+disjunction f a b = case (a, b) of
+  (Nothing, Nothing) -> Nothing
+  (Just x,  Nothing) -> Just x
+  (Nothing, Just y)  -> Just y
+  (Just x, Just y)   -> Just $ f x y
+
+sumMaybe :: Num a => Maybe a -> Maybe a -> Maybe a  
+sumMaybe = disjunction (+)
+
+prodMaybe :: Num a => Maybe a -> Maybe a -> Maybe a  
+prodMaybe = conjunction (*)
+
+
+
+
+
+-- chunk sparse lists according to element index
+
+-- chunkBy q ll = go ll []
+--   where
+--     go _ acc = acc
+--     go (e:es) acc
+--       | q e       = go es (e : acc)
+--       -- | otherwise =
+
+chunks :: Int -> [a] -> [[a]]
+chunks _ [] = []
+chunks n ll = h : chunks n t where
+  (h, t) = splitAt n ll
+
+chunksWhile :: (a -> Bool) -> [a] -> [[a]]
+chunksWhile _ [] = []
+chunksWhile q ll = h : chunksWhile q t where
+  (h, t) = (takeWhile q ll, dropWhile q ll)
+
+chunksWhile' :: (a -> Bool) -> [a] -> [[a]]
+chunksWhile' q = unfoldr genf where
+  genf ll =
+    if null h
+      then Nothing
+      else Just (h, drop (length h) ll)
+    where h = takeWhile q ll 
diff --git a/taco.cabal b/taco.cabal
new file mode 100644
--- /dev/null
+++ b/taco.cabal
@@ -0,0 +1,44 @@
+name:                taco
+version:             0.1.0.0
+synopsis:            Haskell port of the Tensor Algebra COmpiler
+description:         This library provides types and a compiler for tensor expressions.
+homepage:            https://github.com/ocramz/taco-hs#readme
+license:             BSD3
+license-file:        LICENSE
+author:              Marco Zocca
+maintainer:          zocca.marco gmail
+copyright:           2017 Marco Zocca
+category:            Numeric
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+
+library
+  hs-source-dirs:      src
+  ghc-options:         -Wall
+  exposed-modules:     Data.Tensor
+  other-modules:       Data.Shape
+                       Data.Dim
+                       Data.TensorTest
+  build-depends:       base >= 4.7 && < 5
+                     , exceptions
+                     , mtl
+                     , vector
+                     , vector-algorithms
+                     , containers
+  default-language:    Haskell2010
+
+
+
+test-suite test
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Spec.hs
+  build-depends:       base
+                     , taco
+  ghc-options:         -Wall -threaded -rtsopts -with-rtsopts=-N
+  default-language:    Haskell2010
+
+source-repository head
+  type:     git
+  location: https://github.com/ocramz/taco-hs
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,2 @@
+main :: IO ()
+main = putStrLn "Test suite not yet implemented"