diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2012, David Nilsson
+
+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 David Nilsson 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/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/sde-solver.cabal b/sde-solver.cabal
new file mode 100644
--- /dev/null
+++ b/sde-solver.cabal
@@ -0,0 +1,58 @@
+name:                sde-solver
+
+version:             0.1.0.0
+synopsis:            Distributed SDE solver
+description:         
+  This package contains utilities for solving SDE instances in various ways.
+  Basically an 'SDE' instance is solved using some 'SDESolver' working with some distribution mechanism.
+  Results are gathered at the end point of the specified interval.
+  .
+  Included in the package are ways of doing distributed calculations over an MPI cluster,
+  or optionally only using the local solver with built in parallelization.
+  Two SDE instances have been implemented; geometric brownian motion and the Langevin equation,
+  see the haddock documentation for an example.
+  .
+  The main interface is accessible through "Numeric.DSDE" which provides various way of solving generic problems.
+  This module supports either local or distributed calculations in the IO monad and gathering the results as a distribution.
+  Under the surface there is also a working pure implementation for monadic environments, using a pure Mersenne twister PRNG.
+  .
+  Internally there are several abstractions used when dealing with each component building up a solution.
+  Given some 'SDE' and 'SDESolver' instances, it is also required to have some PRNG providing normally distributed numbers.
+  This has been implemented over some specific monads and only results of type 'Double'.
+  All of the internal components are written with polymorphism in mind, acting over some monad instance and generic result types in all cases.
+
+homepage:            https://github.com/davnils/sde-solver
+
+license:             BSD3
+license-file:        LICENSE
+
+author:              David Nilsson
+maintainer:          nilsson.dd+code@gmail.com
+
+category:            Math, Numerical
+build-type:          Simple
+cabal-version:       >=1.8
+
+library
+  exposed-modules:     Numeric.DSDE, Numeric.DSDE.Distribute, Numeric.DSDE.SDESolver, Numeric.DSDE.RNG, Numeric.DSDE.SDE.GeometricBrownian, Numeric.DSDE.SDE.Langevin, Numeric.DSDE.SDE
+  
+  build-depends:
+    base >=4.5 && < 5,
+    cereal >=0.3 && < 1,
+    cereal-vector >=0.2 && < 1,
+    ghc-prim >=0.2 && < 1,
+    haskell-mpi >=1.2 && < 2,
+    mersenne-random-pure64 >= 0.2.0.3 && < 1,
+    mtl >=2.1 && < 3,
+    mwc-random >=0.12 && < 1,
+    normaldistribution >= 1.1.0.3 && < 2,
+    parallel >= 3.2 && < 4,
+    vector >=0.10 && < 1
+  
+  hs-source-dirs:      src
+  ghc-options:         -O2
+
+
+source-repository head
+  type: git
+  location: git://github.com/davnils/sde-solver
diff --git a/src/Numeric/DSDE.hs b/src/Numeric/DSDE.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE.hs
@@ -0,0 +1,83 @@
+{-# LANGUAGE ConstraintKinds #-}
+--------------------------------------------------------------------
+-- |
+-- Module : Numeric.DSDE
+--
+-- Main entry point in the distributed SDE solver.
+--
+-- This module provides wrappers for the most common use cases, 
+-- and the user only needs to choose if to perform a distributed calculation
+-- and supply the problem parameters.
+--
+-- Here follows an example of solving the Langevin equation over
+-- the interval [0,2] with the result being written to a file on the main MPI node.
+-- This produces a file with 1000 entries, corresponding to the final value
+-- in all realisations.
+--
+-- @
+--    let r = 1.0
+--        sigma = 0.1
+--        y_0 = 3.0
+--        end = 0.2
+--        stepSize = 0.01
+--        runs = 1000
+--    withSolver Milstein (distribute (Langevin r sigma) y_0 end stepSize runs) >>=
+--      writeResult "out"
+-- @
+
+module Numeric.DSDE (
+  -- * Abstract solver interface
+  distribute,
+  local,
+  SDEResult(..),
+
+  -- ** Supplied SDE Solvers 
+  Milstein(..),
+  EulerMaruyama(..),
+
+  -- * Helper functions
+  withSolver, writeResult
+) where
+
+import qualified Data.Vector.Unboxed as V
+import Numeric.DSDE.Distribute
+import GHC.Conc
+import Numeric.DSDE.RNG
+import Numeric.DSDE.SDE
+import Numeric.DSDE.SDESolver
+import System.IO
+
+-- | distribute: Perform a distributed MPI calculation over a set of hosts supplied through the MPI environment.
+--   Each host is fully utilized by running a number of GHC threads equivalent to the number cores available.
+--   Internally it uses the highly performant MWC PRNG over the IO monad.
+--
+--   local: Skips the MPI runtime and only runs locally on GHC threads.
+distribute, local :: (SDE sde, SDESolver solver)
+  => sde Double    -- ^ SDE Instance
+  -> Double        -- ^ Initial y_0 value
+  -> Double        -- ^ End of interval, evaluating over [0, end]
+  -> Double        -- ^ Step size used by solver
+  -> Int           -- ^ Number of realisations
+  -> solver        -- ^ Solving method
+  -> IO SDEResult  -- ^ Result 
+
+distribute = evalWrapper [Distr MPI]
+local = evalWrapper [] 
+
+-- | Wrapper function providing access to the internal API.
+evalWrapper :: (SDE sde, SDESolver solver) => [DistributeInstance IO] -> sde Double -> Double -> Double -> Double -> Int -> solver -> IO SDEResult
+evalWrapper distr sde y_0 end step runs solver = do
+  let rng = initialize
+  cores <- getNumCapabilities
+  evaluate (distr, Local cores) (sde, solver, rng, IP (End end) y_0 step runs)
+
+-- | Use the specified solver with an IO action.
+withSolver :: (SDESolver solver, Monad m) => solver -> (solver -> m a) ->  m a
+withSolver = flip id
+
+-- | Write a 'SDEResult' to the supplied filename, output as ASCII text.
+writeResult :: FilePath -> SDEResult -> IO ()
+writeResult file result = withFile file WriteMode $ output result
+  where
+  output (Scalar value count) h = hPutStr h $ unwords [show value, show count]
+  output (Distribution samples) h = V.mapM_ (hPrint h) samples
diff --git a/src/Numeric/DSDE/Distribute.hs b/src/Numeric/DSDE/Distribute.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE/Distribute.hs
@@ -0,0 +1,176 @@
+{-# LANGUAGE MultiParamTypeClasses, ConstraintKinds, BangPatterns,
+             ExistentialQuantification, FlexibleInstances, DeriveGeneric,
+             DefaultSignatures, ScopedTypeVariables, TupleSections,
+             FlexibleContexts #-}
+
+module Numeric.DSDE.Distribute where
+
+import Control.Applicative ((<$>))
+import Control.Concurrent
+import Control.Concurrent.MVar
+import Control.Monad.Identity hiding (mapM)
+import Control.Monad.State
+import Control.Parallel
+import Control.Parallel.MPI.Simple 
+import Data.Foldable (fold, foldl')
+import Data.Monoid
+import Data.Serialize (Serialize(..))
+import qualified Data.Vector.Unboxed as V
+import Data.Vector.Serialize
+import GHC.Generics (Generic)
+import Numeric.DSDE.RNG
+import Numeric.DSDE.SDE
+import Numeric.DSDE.SDESolver
+import Prelude hiding (sum, init, map)
+
+-- | Wrapper used by all distributors supplied to the 'evaluate' function.
+data DistributeInstance m = forall a. Distribute a m => Distr a
+
+-- | Container describing the result produced by an SDE solution.
+data SDEResult = Scalar !Double !Int             -- ^ Average of all end-point values, with the number of samples recorded.
+               | Distribution !(V.Vector Double) -- ^ All end-point samples stored as an unboxed vector.
+  deriving (Generic, Show)
+
+-- | Monoid instance used when folding results from multiple sources.
+instance Monoid SDEResult where
+  mempty = Scalar 0 0
+  Scalar a n         `mappend` Scalar b n'        =
+    Scalar (elemSum / fromIntegral s) s
+    where
+    s = n + n'
+    elemSum = a * fromIntegral n + b * fromIntegral n'
+
+  Distribution v     `mappend` Distribution v'    = Distribution $ v V.++ v'
+  Scalar _ _         `mappend` d@(Distribution _) = d
+  d@(Distribution _) `mappend` Scalar _ _         = d
+
+-- | Serialize instance used by MPI.
+instance Serialize SDEResult
+
+-- | Internal abstraction over the choice of specifying either the interval length or the number of steps.
+data Accuracy = End Double
+              | Steps Int
+  deriving (Generic, Show)
+
+-- | Serialize instance used by MPI.
+instance Serialize Accuracy
+
+-- | Set of parameters supplied to solve an SDE problem.
+data InstanceParams = IP {
+  accuracy :: !Accuracy,
+  start :: !Double,
+  deltat :: !Double,
+  simulations :: !Int }
+  deriving (Generic, Show)
+
+-- | Serialize instance used by MPI.
+instance Serialize InstanceParams
+
+-- | MPI cluster distributor.
+data MPI = MPI
+
+-- | Local evaluation using GHC threads
+data Local = Local Int
+
+type SDEConstraint b c g m p = (SDE b, SDESolver c, Parameter p, RNGGen g m p)
+type SDEInstance b c g m p = (b p, c, Maybe Int -> m g, InstanceParams)
+
+-- | Type class indicating the ability to distribute data in some way.
+--   Several distributors may be chained.
+class Monad m => Distribute a m where
+  -- | Inject an SDE instance into the context.
+  inject :: SDEConstraint b c g m p =>
+    a -> SDEInstance b c g m p -> m (SDEInstance b c g m p)
+  -- | Remove an SDE result from the context.
+  remove :: a -> SDEResult -> m SDEResult
+
+-- | Type class indicating ability to solve an SDE problem.
+class Execute a m p where
+  execute :: SDEConstraint b c g m p =>  a -> SDEInstance b c g m p -> m SDEResult
+
+-- | Type class indicating ability to perform a set of actions in an efficient way.
+class Mappable m p where
+  map' :: RNGGen g m p => (Int, Maybe Int -> m g) -> (g -> m b) -> [a] -> m [b] 
+
+-- | Mappable instance for the IO monad. Work is divided using forkIO.
+instance Mappable IO Double where
+  map' (seed, rng) f l = do
+    rand <- rng (Just seed)
+    seeds <- mapM (\_ -> (Just . round <$> getRand rand) >>= rng) l
+    mapM splitWork seeds >>= mapM takeMVar
+    where
+    splitWork rand = do
+      var <- newEmptyMVar
+      forkIO $ f rand >>= putMVar var
+      return var
+
+-- | Mappable instance for the pure State monad, uses 'par' annotations.
+--   This does not perform well in general and needs to be optimized to
+--   compete with the monadic IO instance.
+instance RealFrac a => Mappable (State s) a where
+  map' (seed, rng) f l = (go f seed l >>= sequence)
+    where
+    go f _ [] = return []
+    go f s (_:t) = do
+      worker <- rng (Just s)
+      s' <- round <$> getRand worker
+      rest <- go f s' t
+      return $ f worker `par` f worker : rest
+
+-- | Distribute instance over MPI which defines data transportation.
+instance Distribute MPI IO where
+  inject _ (sde, solver, rng, params) = do
+    init
+    size <- commSize commWorld
+    rank <- commRank commWorld
+    (sde, solver, rng,) <$> case rank of
+      0 -> do
+        let slaveSize = ceiling $ (fromIntegral $ simulations params :: Double) /
+                                   fromIntegral size
+        let slave = params { simulations = slaveSize}
+        bcastSend commWorld 0 slave
+        return slave
+
+      _ -> bcastRecv commWorld 0
+
+  remove _ localResult = do
+    result <- commRank commWorld >>= retrieve
+    finalize
+    return result
+    where
+    retrieve 0 = do
+      clusterResult <- gatherRecv commWorld 0 localResult
+      return $ fold clusterResult
+    retrieve _ =
+      gatherSend commWorld 0 localResult >> return localResult
+
+-- | Generic execute instances over any mappable monad 'm'.
+instance Mappable m Double => Execute Local m Double where
+  execute (Local cores) (!sde, !solver, !rng, params) = do
+    seedRNG <- round <$> (rng Nothing >>= getRand)
+    fold <$> map' (seedRNG, rng) runThread [1..cores]
+    where
+    perThread = ceiling $ (fromIntegral $ simulations params :: Double) /
+                           fromIntegral cores :: Int
+    steps = case accuracy params of
+      End endTime -> floor $ endTime / deltat params
+      Steps n -> n
+
+    runThread rng = Distribution <$> thread rng
+    thread rand = V.mapM (const $ threadEvaluation rand) $ V.replicate perThread (0.0 :: Double)
+    threadEvaluation rand = foldM' (eval rand) (start params) [1..steps]
+    eval rand w_i step = w_iplus1 solver sde rand (fromIntegral step * deltat params) w_i (deltat params)
+
+-- | Evaluate the SDE using the supplied distributors and execution method.
+evaluate :: (Monad m, SDEConstraint b c g m p, Execute e m p) =>
+  ([DistributeInstance m], e) -> SDEInstance b c g m p -> m SDEResult
+evaluate ([], method) input = execute method input
+evaluate (Distr method : other, final) input =
+  inject method input >>= evaluate (other, final) >>= remove method
+
+-- | Monadic strict fold.
+foldM' :: Monad m => (a -> b -> m a) -> a -> [b] -> m a
+foldM' _ z [] = return z
+foldM' func z (x:xs) = do
+  z' <- func z x
+  z' `seq` foldM' func z' xs
diff --git a/src/Numeric/DSDE/RNG.hs b/src/Numeric/DSDE/RNG.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE/RNG.hs
@@ -0,0 +1,51 @@
+{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances,
+             TypeFamilies, IncoherentInstances, ConstraintKinds,
+             FunctionalDependencies, BangPatterns #-}
+
+module Numeric.DSDE.RNG where
+
+import Control.Applicative
+import Control.Monad.Identity
+import Control.Monad.ST
+import Control.Monad.State
+import Data.Maybe
+import qualified Data.Random.Normal as NORMAL
+import qualified Data.Vector.Unboxed as V
+import qualified System.Random.Mersenne.Pure64 as MT
+import Numeric.DSDE.SDE (Parameter)
+import qualified System.Random.MWC as MWC
+import qualified System.Random.MWC.Distributions as MD
+
+-- | Default seed used by instances when the 'initialize' method is provided with 'Nothing'.
+defaultSeed :: Int
+defaultSeed = 0
+
+-- | Typeclass describing a PRNG working in some monad 'm' and generating values of type 'p'.
+class (Monad m, Functor m, Parameter p) => RNGGen g m p | g m -> p where
+  -- | Generate a N(0,1) random number using the supplied generator.
+  getRand :: g -> m p
+  -- | Initialize a generator using some seed or defaulting to a constant.
+  initialize :: Maybe Int -> m g
+
+-- | 'RNGGen' instance for the ST-monadic MWC RNG.
+instance RNGGen (MWC.Gen s) (ST s) Double where
+  getRand = MD.normal 0 1
+  initialize s = MWC.initialize . V.singleton . fromIntegral $ fromMaybe defaultSeed s
+
+-- | 'RNGGen' instance for the IO-monadic MWC RNG.
+instance (MWC.GenIO ~ d) => RNGGen d IO Double where
+  {-# INLINE getRand #-}
+  getRand = MD.normal 0 1
+  {-# INLINE initialize  #-}
+  initialize (Just n) = MWC.initialize . V.singleton . fromIntegral $ n
+  initialize Nothing = MWC.withSystemRandom . MWC.asGenIO $ return
+
+-- | 'RNGGen' instance for the pure monadic Mersenne Twister RNG, operating in monad 'State'.
+instance RNGGen MT.PureMT (State MT.PureMT) Double where
+  {-# INLINE getRand #-}
+  getRand _ = do
+    !(sample, gen') <- NORMAL.normal <$> get
+    put gen'
+    return sample
+  {-# INLINE initialize  #-}
+  initialize = return . MT.pureMT . fromIntegral . fromMaybe defaultSeed
diff --git a/src/Numeric/DSDE/SDE.hs b/src/Numeric/DSDE/SDE.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE/SDE.hs
@@ -0,0 +1,11 @@
+{-# LANGUAGE ConstraintKinds #-}
+
+module Numeric.DSDE.SDE where
+
+-- | Type of result and intermediate variables used when solving an SDE problem.
+type Parameter a = (Floating a, Fractional a, Num a)
+
+-- | Describes an SDE problem. The partial derivate is used by advanced solvers
+--   such as the 'Milstein' method.
+class SDE a where
+  f,g,partgoverparty :: Parameter p => a p -> p -> p -> p
diff --git a/src/Numeric/DSDE/SDE/GeometricBrownian.hs b/src/Numeric/DSDE/SDE/GeometricBrownian.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE/SDE/GeometricBrownian.hs
@@ -0,0 +1,16 @@
+{-# Language BangPatterns #-}
+
+module Numeric.DSDE.SDE.GeometricBrownian where
+
+import Numeric.DSDE.SDE
+
+-- | Geometric Brownian motion with rate r and diffusion sigma.
+data GeometricBrownian p = GB !p !p
+
+instance SDE GeometricBrownian where
+  {-# SPECIALIZE INLINE f :: GeometricBrownian Double -> Double -> Double -> Double #-}
+  f !(GB rate _) _ !w_i = rate *  w_i
+  {-# SPECIALIZE INLINE g :: GeometricBrownian Double -> Double -> Double -> Double #-}
+  g !(GB _ sigma) _ !w_i = sigma * w_i
+  {-# SPECIALIZE INLINE partgoverparty :: GeometricBrownian Double -> Double -> Double -> Double #-}
+  partgoverparty !(GB _ sigma) _ _ = sigma
diff --git a/src/Numeric/DSDE/SDE/Langevin.hs b/src/Numeric/DSDE/SDE/Langevin.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE/SDE/Langevin.hs
@@ -0,0 +1,16 @@
+{-# Language BangPatterns #-}
+
+module Numeric.DSDE.SDE.Langevin where
+
+import Numeric.DSDE.SDE
+
+-- | SDE instance of the Langevin equation with rate r and diffusion sigma.
+data Langevin p = Langevin !p !p
+
+instance SDE Langevin where
+  {-# SPECIALIZE INLINE f :: Langevin Double -> Double -> Double -> Double #-}
+  f !(Langevin r _) _ !w_i = -r * w_i
+  {-# SPECIALIZE INLINE g :: Langevin Double -> Double -> Double -> Double #-}
+  g !(Langevin _ sigma) _ _ = sigma
+  {-# SPECIALIZE INLINE partgoverparty :: Langevin Double -> Double -> Double -> Double #-}
+  partgoverparty _ _ _ = 0
diff --git a/src/Numeric/DSDE/SDESolver.hs b/src/Numeric/DSDE/SDESolver.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/DSDE/SDESolver.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances,
+             BangPatterns, ConstraintKinds #-}
+
+module Numeric.DSDE.SDESolver where
+
+import Numeric.DSDE.SDE.GeometricBrownian
+import Numeric.DSDE.RNG
+import Numeric.DSDE.SDE
+import qualified System.Random.MWC as M
+
+-- | The Euler-Maruyama solving method. Order 1/2.
+data EulerMaruyama = EulerMaruyama
+
+-- | The Milstein solving method. Order 1.
+data Milstein = Milstein
+
+-- | Type class describing a method of solving SDE problems.
+--    Defined by the next value produced in a solving sequence.
+class SDESolver a where
+  w_iplus1 :: (Monad m, SDE sde, RNGGen rng m p, Parameter p) =>
+    a -> sde p -> rng -> p -> p -> p -> m p
+  solverName :: a -> String
+
+instance SDESolver EulerMaruyama where
+  {-# INLINE w_iplus1 #-}
+  {-# SPECIALIZE w_iplus1 :: EulerMaruyama -> GeometricBrownian Double -> M.GenIO -> Double -> Double -> Double -> IO Double #-}
+  w_iplus1 _ !sde !rng !t_i !w_i !deltat = getRand rng >>= \rand -> return $
+                            w_i
+                            + f sde t_i w_i * deltat
+                            + g sde t_i w_i * deltaB rand
+    where deltaB r = sqrt deltat * r
+
+  solverName _ = "Euler-Maruyama"
+
+instance SDESolver Milstein where
+  w_iplus1 _ !sde !rng !t_i !w_i !deltat = getRand rng >>= \rand -> return $
+                            w_i
+                            + f sde t_i w_i * deltat
+                            + g' * deltaB rand
+                            + g'/2 * partgoverparty sde t_i w_i  * (deltaB rand^^(2 :: Integer) - deltat)
+    where
+    deltaB r = sqrt deltat * r
+    g' = g sde t_i w_i
+  solverName _ = "Milstein"