diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,20 @@
+Copyright (c) 2015 Kai Zhang
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be included
+in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
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/benchmarks/bench.hs b/benchmarks/bench.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/bench.hs
@@ -0,0 +1,38 @@
+import Control.Monad (replicateM)
+import Criterion.Main
+import qualified Data.Clustering.Hierarchical as C
+import qualified Data.Vector as V
+import System.Random.MWC
+
+import AI.Clustering.Hierarchical
+import AI.Clustering.Hierarchical.Types ((!))
+
+randSample :: IO [V.Vector Double]
+randSample = do
+    g <- create
+    replicateM 2000 $ uniformVector g 5
+
+main :: IO ()
+main = do
+    xs <- randSample
+    let dists = computeDists euclidean $ V.fromList xs
+        fn i j = dists ! (i,j)
+    defaultMain
+      [ bgroup "AI.Clustering.Hierarchical"
+          [ bench "Average Linkage (n = 10)" $
+                whnf (\x -> hclust Average x fn) $! V.enumFromN 0 10
+          , bench "Average Linkage (n = 100)" $
+                whnf (\x -> hclust Average x fn) $! V.enumFromN 0 100
+          , bench "Average Linkage (n = 1000)" $
+                whnf (\x -> hclust Average x fn) $! V.enumFromN 0 1000
+          ]
+
+      , bgroup "Data.Clustering.Hierarchical"
+          [ bench "Average Linkage (n = 10)" $
+                whnf (\x -> C.dendrogram C.UPGMA x euclidean) $! take 10 xs
+          , bench "Average Linkage (n = 100)" $
+                whnf (\x -> C.dendrogram C.UPGMA x euclidean) $! take 100 xs
+          , bench "Average Linkage (n = 1000)" $
+                whnf (\x -> C.dendrogram C.UPGMA x euclidean) $! take 1000 xs
+          ]
+      ]
diff --git a/clustering.cabal b/clustering.cabal
new file mode 100644
--- /dev/null
+++ b/clustering.cabal
@@ -0,0 +1,69 @@
+-- Initial fastcluster.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name:                clustering
+version:             0.1.0
+synopsis:            fast clustering algorithms
+description:         O(N^2) implementations for a wide range of hierarchical
+                     clustering schemes, including complete linkage, single linkage,
+                     average linkage, weighted linkage, and Ward's method.
+license:             MIT
+license-file:        LICENSE
+author:              Kai Zhang
+maintainer:          kai@kzhang.org
+copyright:           (c) 2015 Kai Zhang
+category:            Math
+build-type:          Simple
+-- extra-source-files:  
+cabal-version:       >=1.10
+
+library
+  exposed-modules:     
+    AI.Clustering.Hierarchical
+    AI.Clustering.Hierarchical.Types
+
+  other-modules:       
+    AI.Clustering.Hierarchical.Internal
+
+  build-depends:
+      base >=4.0 && <5.0
+    , vector
+    , containers
+
+  hs-source-dirs:      src
+  default-language:    Haskell2010
+
+test-suite test
+  type: exitcode-stdio-1.0
+  hs-source-dirs: tests
+  main-is: test.hs
+  other-modules:
+    Test.Hierarchical
+
+  default-language:    Haskell2010
+  build-depends: 
+      base
+    , mwc-random
+    , vector
+    , tasty
+    , tasty-hunit
+    , clustering
+    , hierarchical-clustering
+
+benchmark bench
+  type: exitcode-stdio-1.0
+  hs-source-dirs: benchmarks
+  main-is: bench.hs
+
+  default-language:    Haskell2010
+  build-depends: 
+      base
+    , criterion
+    , mwc-random
+    , vector
+    , clustering
+    , hierarchical-clustering
+
+source-repository  head
+  type: git
+  location: https://github.com/kaizhang/clustering.git
diff --git a/src/AI/Clustering/Hierarchical.hs b/src/AI/Clustering/Hierarchical.hs
new file mode 100644
--- /dev/null
+++ b/src/AI/Clustering/Hierarchical.hs
@@ -0,0 +1,64 @@
+{-# LANGUAGE FlexibleContexts #-}
+--------------------------------------------------------------------------------
+-- |
+-- Module      :  $Header$
+-- Description :  <optional short text displayed on contents page>
+-- Copyright   :  (c) Kai Zhang
+-- License     :  MIT
+
+-- Maintainer  :  kai@kzhang.org
+-- Stability   :  experimental
+-- Portability :  portable
+
+-- <module description starting at first column>
+--------------------------------------------------------------------------------
+
+module AI.Clustering.Hierarchical
+    ( Dendrogram(..)
+    , size
+    , cutAt
+    , members
+    , Metric(..)
+    , hclust
+    , computeDists
+    , euclidean
+    ) where
+
+import Control.Applicative ((<$>))
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+import AI.Clustering.Hierarchical.Internal
+import AI.Clustering.Hierarchical.Types
+
+data Metric = Single    -- ^ Single linkage, $d(A,B) = min_{a \in A, b \in B} d(a,b)$.
+            | Complete  -- ^ Complete linkage, $d(A,B) = max_{a \in A, b \in B} d(a,b)$.
+            | Average   -- ^ Average linkage, $d(A,B) = \frac{\sum_{a \in A}\sum_{b \in B}d(a,b)}{|A||B|}$.
+            | Weighted  -- ^ Weighted linkage
+            | Ward      -- ^ Ward's method
+            | Centroid  -- ^ Centroid linkage, not implemented
+            | Median    -- ^ Median linkage, not implemented
+
+hclust :: G.Vector v a => Metric -> v a -> DistFn a -> Dendrogram a
+hclust method xs f = label <$> nnChain dists fn
+  where
+    dists = computeDists f xs
+    label i = xs G.! i
+    fn = case method of
+        Single -> single
+        Complete -> complete
+        Average -> average
+        Weighted -> weighted
+        Ward -> ward
+        _ -> error "Not implemented"
+
+computeDists :: G.Vector v a => DistFn a -> v a -> DistanceMat
+computeDists f vec = DistanceMat n . U.fromList . flip concatMap [0..n-1] $ \i ->
+    flip map [i+1..n-1] $ \j -> f (vec `G.unsafeIndex` i) (vec `G.unsafeIndex` j)
+  where
+    n = G.length vec
+{-# INLINE computeDists #-}
+
+euclidean :: G.Vector v Double => DistFn (v Double)
+euclidean xs ys = sqrt $ G.sum $ G.zipWith (\x y -> (x-y)**2) xs ys
+{-# INLINE euclidean #-}
diff --git a/src/AI/Clustering/Hierarchical/Internal.hs b/src/AI/Clustering/Hierarchical/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/AI/Clustering/Hierarchical/Internal.hs
@@ -0,0 +1,122 @@
+module AI.Clustering.Hierarchical.Internal
+    ( nnChain
+    , single
+    , complete
+    , average
+    , weighted
+    , ward
+    ) where
+
+import Control.Monad (forM_, when)
+import qualified Data.Map as M
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as UM
+
+import AI.Clustering.Hierarchical.Types
+
+type ActiveNodeSet = M.Map Int (Dendrogram Int)
+type DistUpdateFn = Int -> Int -> ActiveNodeSet -> DistanceMat -> DistanceMat
+
+-- | nearest neighbor chain algorithm
+nnChain :: DistanceMat -> DistUpdateFn -> Dendrogram Int
+nnChain (DistanceMat n dist) fn = go (DistanceMat n $ U.force dist) initSet []
+  where
+    go ds activeNodes chain@(b:a:rest)
+        | M.size activeNodes == 1 = head . M.elems $ activeNodes
+        | c == a = go ds' activeNodes' rest
+        | otherwise = go ds activeNodes $ c : chain
+      where
+        (c,d) = nearestNeighbor ds b a activeNodes
+        activeNodes' = M.insert hi (Branch (size1+size2) d c1 c2)
+                     . M.delete lo $ activeNodes
+        ds' = fn lo hi activeNodes ds
+        c1 = M.findWithDefault undefined lo activeNodes
+        c2 = M.findWithDefault undefined hi activeNodes
+        size1 = size c1
+        size2 = size c2
+        (lo,hi) = if a <= b then (a,b) else (b,a)
+    go ds activeNodes _ = go ds activeNodes [b,a]
+      where
+        a = fst $ M.elemAt 0 activeNodes
+        b = fst $ nearestNeighbor ds a (-1) activeNodes
+    initSet = M.fromList . map (\i -> (i, Leaf i)) $ [0..n-1]
+{-# INLINE nnChain #-}
+
+nearestNeighbor :: DistanceMat -> Int -> Int -> M.Map Int (Dendrogram Int) -> (Int, Double)
+nearestNeighbor dist i preference = M.foldlWithKey' f (-1,1/0)
+  where
+    f (x,d) j _ | i == j = (x,d)  -- skip
+                | d' < d = (j,d')
+                | d' == d && j == preference = (j,d')
+                | otherwise = (x,d)
+      where d' = dist ! (i,j)
+{-# INLINE nearestNeighbor #-}
+
+-- | all update functions perform destructive updates, and hence should not be
+-- called outside this module
+
+-- | single linkage update formula
+single :: DistUpdateFn
+single lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do
+    v <- U.unsafeThaw dist
+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do
+        d_lo_i <- UM.unsafeRead v $ idx n i lo
+        d_hi_i <- UM.unsafeRead v $ idx n i hi
+        UM.unsafeWrite v (idx n i hi) $ min d_lo_i d_hi_i
+    return v
+{-# INLINE single #-}
+
+-- | complete linkage update formula
+complete :: DistUpdateFn
+complete lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do
+    v <- U.unsafeThaw dist
+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do
+        d_lo_i <- UM.unsafeRead v $ idx n i lo
+        d_hi_i <- UM.unsafeRead v $ idx n i hi
+        UM.unsafeWrite v (idx n i hi) $ max d_lo_i d_hi_i
+    return v
+{-# INLINE complete #-}
+
+-- | average linkage update formula
+average :: DistUpdateFn
+average lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do
+    v <- U.unsafeThaw dist
+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do
+        d_lo_i <- UM.unsafeRead v $ idx n i lo
+        d_hi_i <- UM.unsafeRead v $ idx n i hi
+        UM.unsafeWrite v (idx n i hi) $ f1 * d_lo_i + f2 * d_hi_i
+    return v
+  where
+    s1 = fromIntegral . size . M.findWithDefault undefined lo $ nodeset
+    s2 = fromIntegral . size . M.findWithDefault undefined hi $ nodeset
+    f1 = s1 / (s1+s2)
+    f2 = s2 / (s1+s2)
+{-# INLINE average #-}
+
+-- | complete linkage update formula
+weighted :: DistUpdateFn
+weighted lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do
+    v <- U.unsafeThaw dist
+    forM_ (M.keys nodeset) $ \i -> when (i/= hi && i/=lo) $ do
+        d_lo_i <- UM.unsafeRead v $ idx n i lo
+        d_hi_i <- UM.unsafeRead v $ idx n i hi
+        UM.unsafeWrite v (idx n i hi) $ (d_lo_i + d_hi_i) / 2
+    return v
+{-# INLINE weighted #-}
+
+-- | ward linkage update formula
+ward :: DistUpdateFn
+ward lo hi nodeset (DistanceMat n dist) = DistanceMat n $ U.create $ do
+    v <- U.unsafeThaw dist
+    d_lo_hi <- UM.unsafeRead v $ idx n lo hi
+    forM_ (M.toList nodeset) $ \(i,t) -> when (i/= hi && i/=lo) $ do
+        let s3 = fromIntegral . size $ t
+        d_lo_i <- UM.unsafeRead v $ idx n i lo
+        d_hi_i <- UM.unsafeRead v $ idx n i hi
+        UM.unsafeWrite v (idx n i hi) $
+            sqrt $ ((s1+s3)*d_lo_i + (s2+s3)*d_hi_i - s3*d_lo_hi) / (s1+s2+s3)
+    return v
+  where
+    s1 = fromIntegral . size . M.findWithDefault undefined lo $ nodeset
+    s2 = fromIntegral . size . M.findWithDefault undefined hi $ nodeset
+{-# INLINE ward #-}
diff --git a/src/AI/Clustering/Hierarchical/Types.hs b/src/AI/Clustering/Hierarchical/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/AI/Clustering/Hierarchical/Types.hs
@@ -0,0 +1,55 @@
+module AI.Clustering.Hierarchical.Types
+    ( Distance
+    , DistFn
+    , Size
+    , Dendrogram(..)
+    , size
+    , cutAt
+    , members
+    , DistanceMat(..)
+    , (!)
+    , idx
+    ) where
+
+import Data.Bits (shiftR)
+import qualified Data.Vector.Unboxed as U
+
+type Distance = Double
+type DistFn a = a -> a -> Distance
+type Size = Int
+
+data Dendrogram a = Leaf !a
+                  | Branch !Size !Distance !(Dendrogram a) !(Dendrogram a)
+    deriving (Show)
+
+instance Functor Dendrogram where
+    fmap f (Leaf x) = Leaf $ f x
+    fmap f (Branch n d l r) = Branch n d (fmap f l) $ fmap f r
+
+size :: Dendrogram a -> Int
+size (Leaf _) = 1
+size (Branch n _ _ _) = n
+{-# INLINE size #-}
+
+cutAt :: Dendrogram a -> Distance -> [Dendrogram a]
+cutAt dendro th = go [] dendro
+  where
+    go acc x@(Leaf _) = x : acc
+    go acc x@(Branch _ d l r) | d <= th = x : acc
+                              | otherwise = go (go acc r) l
+
+members :: Dendrogram a -> [a]
+members (Leaf x) = [x]
+members (Branch _ _ l r) = members l ++ members r
+
+-- upper triangular matrix
+data DistanceMat = DistanceMat !Int !(U.Vector Double) deriving (Show)
+
+(!) :: DistanceMat -> (Int, Int) -> Double
+(!) (DistanceMat n v) (i',j') = v U.! idx n i' j'
+{-# INLINE (!) #-}
+
+idx :: Int -> Int -> Int -> Int
+idx n i j | i <= j = (i * (2 * n - i - 3)) `shiftR` 1 + j - 1
+          | otherwise = (j * (2 * n - j - 3)) `shiftR` 1 + i - 1
+{-# INLINE idx #-}
diff --git a/tests/Test/Hierarchical.hs b/tests/Test/Hierarchical.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/Hierarchical.hs
@@ -0,0 +1,64 @@
+module Test.Hierarchical
+    ( tests
+    ) where
+
+import Control.Monad
+import qualified Data.Clustering.Hierarchical as C
+import qualified Data.Vector as V
+import System.Random.MWC
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import AI.Clustering.Hierarchical
+
+tests :: TestTree
+tests = testGroup "Hierarchical:"
+    [ testCase "Single Linkage" testSingle
+    , testCase "Complete Linkage" testComplete
+    , testCase "Average Linkage" testAverage
+    , testCase "Weighted Linkage" testWeighted
+    ]
+
+randSample :: IO [V.Vector Double]
+randSample = do
+    g <- create
+    replicateM 500 $ uniformVector g 5
+
+isEqual :: Eq a => Dendrogram a -> C.Dendrogram a -> Bool
+isEqual (Leaf x) (C.Leaf x') = x == x'
+isEqual (Branch _ d x y) (C.Branch d' x' y') = abs (d - d') < 1e-8 &&
+    ((isEqual x x' && isEqual y y') || (isEqual x y' && isEqual y x'))
+isEqual _ _ = False
+
+testSingle :: Assertion
+testSingle = do
+    xs <- randSample
+    let true = C.dendrogram C.SingleLinkage xs euclidean
+        test = hclust Single (V.fromList xs) euclidean
+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $
+        isEqual test true
+
+testComplete :: Assertion
+testComplete = do
+    xs <- randSample
+    let true = C.dendrogram C.CompleteLinkage xs euclidean
+        test = hclust Complete (V.fromList xs) euclidean
+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $
+        isEqual test true
+
+testAverage :: Assertion
+testAverage = do
+    xs <- randSample
+    let true = C.dendrogram C.UPGMA xs euclidean
+        test = hclust Average (V.fromList xs) euclidean
+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $
+        isEqual test true
+
+testWeighted :: Assertion
+testWeighted = do
+    xs <- randSample
+    let true = C.dendrogram C.FakeAverageLinkage xs euclidean
+        test = hclust Weighted (V.fromList xs) euclidean
+    assertBool (unlines ["Expect: ", show true, "But see: ", show test]) $
+        isEqual test true
+
diff --git a/tests/test.hs b/tests/test.hs
new file mode 100644
--- /dev/null
+++ b/tests/test.hs
@@ -0,0 +1,7 @@
+import Test.Tasty
+
+import qualified Test.Hierarchical as Hierarchical
+
+main :: IO ()
+main = defaultMain $ testGroup "Main"
+    [ Hierarchical.tests ]