80 files changed, 2542 insertions, 68 deletions

diff --git a/.circleci/config.yml b/.circleci/config.yml
new file mode 100644
index 0000000..68fa538
--- /dev/null
+++ b/.circleci/config.yml
@@ -0,0 +1,33 @@
+version: 2
+jobs:
+  build:
+    docker:
+      - image: terrorjack/vanilla:haskell
+    steps:
+      - checkout
+      - restore_cache:
+          key: stack-deps-eek-{{ checksum "stack.yaml" }}
+      - run:
+          name: Update packages
+          command: nix-channel --update
+      - run:
+          name: Building dependencies
+          command: |
+            nix-env -f nixpkgs.nix -iA stack
+            stack --nix --no-terminal --nix build --only-snapshot --prefetch --no-haddock --test --bench --jobs=1
+      - save_cache:
+          paths:
+            - "~/.stack"
+          key: stack-deps-eek-{{ checksum "stack.yaml" }}
+      - run:
+          name: Building
+          command: |
+            stack --nix --no-terminal --nix build --pedantic
+      - run:
+          name: Building tests
+          command: |
+            stack --nix test --pedantic --no-run-tests
+      - run:
+          name: Running tests
+          command: |
+            stack --nix test
diff --git a/INSTALL.md b/INSTALL.md
index 9f5761f..99ddc38 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -1,6 +1,15 @@
 # hmatrix installation
 
-This package requires [GHC 7.8](http://www.haskell.org/ghc), [cabal-install](http://www.haskell.org/haskellwiki/Cabal-Install) (available in the [Haskell Platform](http://hackage.haskell.org/platform)), and the development packages for BLAS/[LAPACK](http://www.netlib.org/lapack) and [GSL](http://www.gnu.org/software/gsl).
+*WARNING*: these instructions are now very out of date. However, there
+is a CI script which shows what does work; you are advised to look at
+that and emulate its installation steps.
+
+This package requires [GHC 7.8](http://www.haskell.org/ghc),
+[cabal-install](http://www.haskell.org/haskellwiki/Cabal-Install)
+(available in the [Haskell
+Platform](http://hackage.haskell.org/platform)), and the development
+packages for BLAS/[LAPACK](http://www.netlib.org/lapack) and
+[GSL](http://www.gnu.org/software/gsl).
 
 ## Linux ##################################################
 
diff --git a/README.md b/README.md
index cdbcf85..de7f523 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,22 @@
 A Haskell library for numerical computation
 -------------------------------------------
 
-A purely functional interface to linear algebra and other numerical algorithms, internally implemented using [LAPACK][lapack], [BLAS][blas], and [GSL][gsl].
+A purely functional interface to linear algebra and other numerical
+algorithms, internally implemented using [LAPACK][lapack],
+[BLAS][blas], [GSL][gsl] and [SUNDIALS][sundials].
 
-This package includes matrix decompositions (eigensystems, singular values, Cholesky, QR, etc.), linear solvers, numeric integration, root finding, etc.
+This package includes matrix decompositions (eigensystems, singular
+values, Cholesky, QR, etc.), linear solvers, numeric integration, root
+finding, etc.
 
-- [What's new][changes] in version 0.17 (July 2015)
+- [What's new][changes] in version 0.19 (April 2018). This is not
+  intended to be a breaking change but a lot of modules have been
+  modified to ensure that continuous integration is green.
+
+- Support for SUNDIALS has been added. It should be possible to
+  replace Numeric.GSL.ODE with Numeric.Sundials.ARKode.ODE and have
+  your program work as before bearing in mind that the methods and
+  error control might differ (even for those with the same names!).
 
 - [Code examples][examples]
 
@@ -26,6 +37,7 @@ Contributions, suggestions, and bug reports are welcome!
 [lapack]: http://www.netlib.org/lapack/
 [blas]: http://www.netlib.org/blas/
 [gsl]: http://www.gnu.org/software/gsl/
+[sundials]: https://computation.llnl.gov/projects/sundials
 
 [tutorial]: http://dis.um.es/profesores/alberto/material/hmatrix.pdf
 [installation]: https://github.com/AlbertoRuiz/hmatrix/blob/master/INSTALL.md
diff --git a/examples/ButcherTableau.hs b/examples/ButcherTableau.hs
new file mode 100644
index 0000000..e667ef1
--- /dev/null
+++ b/examples/ButcherTableau.hs
@@ -0,0 +1,47 @@
+{-# OPTIONS_GHC -Wall #-}
+
+import           Numeric.Sundials.ARKode.ODE
+import           Numeric.LinearAlgebra
+
+import           Data.List (intercalate)
+
+import           Text.PrettyPrint.HughesPJClass
+
+
+butcherTableauTex :: ButcherTable -> String
+butcherTableauTex (ButcherTable m c b b2) =
+  render $
+  vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}")
+       , us
+       , text "\\hline"
+       , text bs <+> text "\\\\"
+       , text b2s <+> text "\\\\"
+       , text "\\end{array}"
+       ]
+  where
+    n = rows m
+    rs = toLists m
+    ss = map (\r -> intercalate " & " $ map show r) rs
+    ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss
+    us = vcat $ map (\r -> text r <+> text "\\\\") ts
+    bs  = " & " ++ (intercalate " & " $ map show $ toList b)
+    b2s = " & " ++ (intercalate " & " $ map show $ toList b2)
+
+main :: IO ()
+main = do
+
+  let res = butcherTable (SDIRK_2_1_2 undefined)
+  putStrLn $ show res
+  putStrLn $ butcherTableauTex res
+
+  let resA = butcherTable (KVAERNO_4_2_3 undefined)
+  putStrLn $ show resA
+  putStrLn $ butcherTableauTex resA
+
+  let resB = butcherTable (SDIRK_5_3_4 undefined)
+  putStrLn $ show resB
+  putStrLn $ butcherTableauTex resB
+
+  let resC = butcherTable (FEHLBERG_6_4_5 undefined)
+  putStrLn $ show resC
+  putStrLn $ butcherTableauTex resC
diff --git a/examples/ChangeLog.md b/examples/ChangeLog.md
new file mode 100644
index 0000000..a7e4c4c
--- /dev/null
+++ b/examples/ChangeLog.md
@@ -0,0 +1,5 @@
+# Revision history for hmatrix-examples
+
+## 0.1.0.0  -- YYYY-mm-dd
+
+* First version. Released on an unsuspecting world.
diff --git a/examples/LICENSE b/examples/LICENSE
new file mode 100644
index 0000000..18e5956
--- /dev/null
+++ b/examples/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2018, Dominic Steinitz
+
+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 Dominic Steinitz 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/examples/Setup.hs b/examples/Setup.hs
new file mode 100644
index 0000000..9a994af
--- /dev/null
+++ b/examples/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/examples/examples.cabal b/examples/examples.cabal
new file mode 100644
index 0000000..07ec571
--- /dev/null
+++ b/examples/examples.cabal
@@ -0,0 +1,29 @@
+name:                examples
+version:             0.19.0.0
+synopsis:            Example usage of the various hmatrix packages
+homepage:            https://github.com/albertoruiz/hmatrix
+license:             BSD3
+license-file:        LICENSE
+author:              Dominic Steinitz
+maintainer:          dominic@steinitz.org
+copyright:           Dominic Steinitz, Novadiscovery
+category:            Math
+build-type:          Simple
+extra-source-files:  ChangeLog.md
+cabal-version:       >=1.10
+
+executable sundials
+  main-is:             sundials.hs
+  build-depends:       base >=4.10 && <4.11,
+                       hmatrix,
+                       hmatrix-sundials,
+                       hmatrix-gsl
+  default-language:    Haskell2010
+
+executable butcherTableau
+  main-is:             ButcherTableau.hs
+  build-depends:       base >=4.10 && <4.11,
+                       hmatrix,
+                       hmatrix-sundials,
+                       pretty
+  default-language:    Haskell2010
diff --git a/examples/sundials.hs b/examples/sundials.hs
new file mode 100644
index 0000000..99f662d
--- /dev/null
+++ b/examples/sundials.hs
@@ -0,0 +1,56 @@
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
+{-# LANGUAGE ViewPatterns #-}
+
+import           Numeric.Sundials.ARKode.ODE
+import           Numeric.LinearAlgebra
+import           Graphics.Plot
+
+vanderpol mu = do
+    let xdot nu _t [x,v] = [v, -x + nu * v * (1-x*x)]
+        xdot _ _ _ = error "vanderpol RHS not defined"
+        ts = linspace 1000 (0,50)
+        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
+    mplot (ts : sol)
+    mplot sol
+
+
+harmonic w d = do
+    let xdot u dd _t [x,v] = [v, a*x + b*v] where a = -u*u; b = -2*dd*u
+        xdot _ _ _ _ = error "harmonic RHS not defined"
+        ts = linspace 100 (0,20)
+        sol = odeSolve (xdot w d) [1,0] ts
+    mplot (ts : toColumns sol)
+
+
+kepler v a = mplot (take 2 $ toColumns sol) where
+    xdot _t [x,y,vx,vy] = [vx,vy,x*k,y*k]
+        where g=1
+              k=(-g)*(x*x+y*y)**(-1.5)
+    xdot _ _ = error "kepler RHS not defined"
+    ts = linspace 100 (0,30)
+    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
+    degree = pi/180
+
+
+main = do
+    vanderpol 2
+    harmonic 1 0
+    harmonic 1 0.1
+    kepler 0.3 60
+    kepler 0.4 70
+    vanderpol' 2
+
+-- example of odeSolveV with jacobian
+vanderpol' mu = do
+    let xdot nu _t (toList->[x,v]) = fromList [v, -x + nu * v * (1-x*x)]
+        xdot _ _ _ = error "vanderpol' RHS not defined"
+        jac _ (toList->[x,v]) = (2><2) [ 0          ,          1
+                                       , -1-2*x*v*mu, mu*(1-x**2) ]
+        jac _ _ = error "vanderpol' Jacobian not defined"
+        ts = linspace 1000 (0,50)
+        hi = pure $ (ts!1 - ts!0) / 100.0
+        sol = toColumns $ odeSolveV (SDIRK_5_3_4 jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
+    mplot sol
+
+
diff --git a/nixpkgs.nix b/nixpkgs.nix
new file mode 100644
index 0000000..77c89c6
--- /dev/null
+++ b/nixpkgs.nix
@@ -0,0 +1 @@
+import (fetchTarball "https://github.com/nixos/nixpkgs/archive/f93d4b6181f5c28b20d63a4d76e182511369c1bf.tar.gz")
diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG
index 11a57dd..fd1c171 100644
--- a/packages/base/CHANGELOG
+++ b/packages/base/CHANGELOG
@@ -2,15 +2,15 @@
 --------
 
     * Many new functions and instances in the Static module
-    
+
     * meanCov and gaussianSample use Herm type
 
     * thinQR, thinRQ
-    
+
     * compactSVDTol
 
     * unitary changed to normalize, also admits Vector (Complex Double)
-  
+
 0.17.0.0
 --------
 
@@ -288,4 +288,3 @@
     * added NFData instances for Matrix and Vector.
 
     * liftVector, liftVector2 replaced by mapVector, zipVector.
-
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 6368683..1380a36 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.18.2.0
+Version:            0.19.0.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index 5fe7796..6027c46 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -4,6 +4,8 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE TypeFamilies #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 {- |
 Module      :  Internal.Algorithms
diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs
index cc10ad8..29edd35 100644
--- a/packages/base/src/Internal/CG.hs
+++ b/packages/base/src/Internal/CG.hs
@@ -1,6 +1,8 @@
 {-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
 {-# LANGUAGE RecordWildCards #-}
 
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
 module Internal.CG(
     cgSolve, cgSolve',
     CGState(..), R, V
diff --git a/packages/base/src/Internal/Chain.hs b/packages/base/src/Internal/Chain.hs
index f87eb02..4000c2b 100644
--- a/packages/base/src/Internal/Chain.hs
+++ b/packages/base/src/Internal/Chain.hs
@@ -1,5 +1,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Internal.Chain
diff --git a/packages/base/src/Internal/Container.hs b/packages/base/src/Internal/Container.hs
index cdcdad0..41b8214 100644
--- a/packages/base/src/Internal/Container.hs
+++ b/packages/base/src/Internal/Container.hs
@@ -5,6 +5,8 @@
 {-# LANGUAGE FunctionalDependencies #-}
 {-# LANGUAGE UndecidableInstances #-}
 
+{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Internal.Container
diff --git a/packages/base/src/Internal/Devel.hs b/packages/base/src/Internal/Devel.hs
index 3887663..f72d8aa 100644
--- a/packages/base/src/Internal/Devel.hs
+++ b/packages/base/src/Internal/Devel.hs
@@ -54,6 +54,7 @@ check msg f = do
 
 -- | postfix error code check
 infixl 0 #|
+(#|) :: IO CInt -> String -> IO ()
 (#|) = flip check
 
 -- | Error capture and conversion to Maybe
diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs
index eb3a25b..2e330ee 100644
--- a/packages/base/src/Internal/Element.hs
+++ b/packages/base/src/Internal/Element.hs
@@ -4,6 +4,8 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Matrix
@@ -31,6 +33,7 @@ import Data.List.Split(chunksOf)
 import Foreign.Storable(Storable)
 import System.IO.Unsafe(unsafePerformIO)
 import Control.Monad(liftM)
+import Foreign.C.Types(CInt)
 
 -------------------------------------------------------------------
 
@@ -53,8 +56,10 @@ instance (Show a, Element a) => (Show (Matrix a)) where
     show m | rows m == 0 || cols m == 0 = sizes m ++" []"
     show m = (sizes m++) . dsp . map (map show) . toLists $ m
 
+sizes :: Matrix t -> [Char]
 sizes m = "("++show (rows m)++"><"++show (cols m)++")\n"
 
+dsp :: [[[Char]]] -> [Char]
 dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
     where
         mt = transpose as
@@ -73,6 +78,7 @@ instance (Element a, Read a) => Read (Matrix a) where
               rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
 
 
+breakAt :: Eq a => a -> [a] -> ([a], [a])
 breakAt c l = (a++[c],tail b) where
     (a,b) = break (==c) l
 
@@ -88,7 +94,8 @@ data Extractor
     | Drop Int
     | DropLast Int
   deriving Show
-  
+
+ppext :: Extractor -> [Char]
 ppext All = ":"
 ppext (Range a 1 c) = printf "%d:%d" a c
 ppext (Range a b c) = printf "%d:%d:%d" a b c
@@ -128,10 +135,14 @@ ppext (DropLast n) = printf "DropLast %d" n
 infixl 9 ??
 (??)  :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t
 
+minEl :: Vector CInt -> CInt
 minEl = toScalarI Min
+maxEl :: Vector CInt -> CInt
 maxEl = toScalarI Max
+cmodi :: Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt
 cmodi = vectorMapValI ModVS
 
+extractError :: Matrix t1 -> (Extractor, Extractor) -> t
 extractError m (e1,e2)= error $ printf "can't extract (%s,%s) from matrix %dx%d" (ppext e1::String) (ppext e2::String) (rows m) (cols m)
 
 m ?? (Range a s b,e) | s /= 1 = m ?? (Pos (idxs [a,a+s .. b]), e)
@@ -232,8 +243,10 @@ disp = putStr . dispf 2
 fromBlocks :: Element t => [[Matrix t]] -> Matrix t
 fromBlocks = fromBlocksRaw . adaptBlocks
 
+fromBlocksRaw :: Element t => [[Matrix t]] -> Matrix t
 fromBlocksRaw mms = joinVert . map joinHoriz $ mms
 
+adaptBlocks :: Element t => [[Matrix t]] -> [[Matrix t]]
 adaptBlocks ms = ms' where
     bc = case common length ms of
           Just c -> c
@@ -486,6 +499,9 @@ liftMatrix2Auto f m1 m2
     m2' = conformMTo (r,c) m2
 
 -- FIXME do not flatten if equal order
+lM :: (Storable t, Element t1, Element t2)
+   => (Vector t1 -> Vector t2 -> Vector t)
+   -> Matrix t1 -> Matrix t2 -> Matrix t
 lM f m1 m2 = matrixFromVector
                 RowMajor
                 (max' (rows m1) (rows m2))
@@ -504,6 +520,7 @@ compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
 
 ------------------------------------------------------------
 
+toBlockRows :: Element t => [Int] -> Matrix t -> [Matrix t]
 toBlockRows [r] m
     | r == rows m = [m]
 toBlockRows rs m
@@ -513,6 +530,7 @@ toBlockRows rs m
     szs = map (* cols m) rs
     g k = (k><0)[]
 
+toBlockCols :: Element t => [Int] -> Matrix t -> [Matrix t]
 toBlockCols [c] m | c == cols m = [m]
 toBlockCols cs m = map trans . toBlockRows cs . trans $ m
 
@@ -576,7 +594,7 @@ Just (3><3)
 mapMatrixWithIndexM
   :: (Element a, Storable b, Monad m) =>
       ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
-mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m 
+mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m
     where
       c = cols m
 
@@ -598,4 +616,3 @@ mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
 
 mapMatrix :: (Element a, Element b) => (a -> b) -> Matrix a -> Matrix b
 mapMatrix f = liftMatrix (mapVector f)
-
diff --git a/packages/base/src/Internal/IO.hs b/packages/base/src/Internal/IO.hs
index a899cfd..b0f5606 100644
--- a/packages/base/src/Internal/IO.hs
+++ b/packages/base/src/Internal/IO.hs
@@ -20,7 +20,7 @@ import Internal.Devel
 import Internal.Vector
 import Internal.Matrix
 import Internal.Vectorized
-import Text.Printf(printf)
+import Text.Printf(printf, PrintfArg, PrintfType)
 import Data.List(intersperse,transpose)
 import Data.Complex
 
@@ -78,12 +78,18 @@ disps d x = sdims x ++ "  " ++ formatScaled d x
 dispf :: Int -> Matrix Double -> String
 dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
 
+sdims :: Matrix t -> [Char]
 sdims x = show (rows x) ++ "x" ++ show (cols x)
 
+formatFixed :: (Show a, Text.Printf.PrintfArg t, Element t)
+            => a -> Matrix t -> String
 formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x
 
+isInt :: Matrix Double -> Bool
 isInt = all lookslikeInt . toList . flatten
 
+formatScaled :: (Text.Printf.PrintfArg b, RealFrac b, Floating b, Num t, Element b, Show t)
+             => t -> Matrix b -> [Char]
 formatScaled dec t = "E"++show o++"\n" ++ ss
     where ss = format " " (printf fmt. g) t
           g x | o >= 0    = x/10^(o::Int)
@@ -133,14 +139,18 @@ showComplex d (a:+b)
     s2 = if b<0 then "-" else ""
     s3 = if b<0 then "-" else "+"
 
+shcr :: (Show a, Show t1, Text.Printf.PrintfType t, Text.Printf.PrintfArg t1, RealFrac t1)
+     => a -> t1 -> t
 shcr d a | lookslikeInt a = printf "%.0f" a
          | otherwise      = printf ("%."++show d++"f") a
 
-
+lookslikeInt :: (Show a, RealFrac a) => a -> Bool
 lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
    where shx = show x
 
+isZero :: Show a => a -> Bool
 isZero x = show x `elem` ["0.0","-0.0"]
+isOne :: Show a => a -> Bool
 isOne  x = show x `elem` ["1.0","-1.0"]
 
 -- | Pretty print a complex matrix with at most n decimal digits.
@@ -168,6 +178,6 @@ loadMatrix f = do
       else
         return (reshape c v)
 
-
+loadMatrix' :: FilePath -> IO (Maybe (Matrix Double))
 loadMatrix' name = mbCatch (loadMatrix name)
 
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index e306454..64cf2f5 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -1,6 +1,8 @@
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE ViewPatterns #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.LinearAlgebra.LAPACK
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index 2856ec2..5436e59 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -57,19 +57,24 @@ cols :: Matrix t -> Int
 cols = icols
 {-# INLINE cols #-}
 
+size :: Matrix t -> (Int, Int)
 size m = (irows m, icols m)
 {-# INLINE size #-}
 
+rowOrder :: Matrix t -> Bool
 rowOrder m = xCol m == 1 || cols m == 1
 {-# INLINE rowOrder #-}
 
+colOrder :: Matrix t -> Bool
 colOrder m = xRow m == 1 || rows m == 1
 {-# INLINE colOrder #-}
 
+is1d :: Matrix t -> Bool
 is1d (size->(r,c)) = r==1 || c==1
 {-# INLINE is1d #-}
 
 -- data is not contiguous
+isSlice :: Storable t => Matrix t -> Bool
 isSlice m@(size->(r,c)) = r*c < dim (xdat m)
 {-# INLINE isSlice #-}
 
@@ -136,16 +141,20 @@ instance Storable t => TransArray (Matrix t)
     {-# INLINE applyRaw #-}
 
 infixr 1 #
+(#) :: TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
 a # b = apply a b
 {-# INLINE (#) #-}
 
+(#!) :: (TransArray c, TransArray c1) => c1 -> c -> Trans c1 (Trans c (IO r)) -> IO r
 a #! b = a # b # id
 {-# INLINE (#!) #-}
 
 --------------------------------------------------------------------------------
 
+copy :: Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
 copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])
 
+extractAll :: Element t => MatrixOrder -> Matrix t -> Matrix t
 extractAll ord m = unsafePerformIO (copy ord m)
 
 {- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
@@ -224,11 +233,13 @@ m@Matrix {irows = r, icols = c} @@> (i,j)
 {-# INLINE (@@>) #-}
 
 --  Unsafe matrix access without range checking
+atM' :: Storable t => Matrix t -> Int -> Int -> t
 atM' m i j = xdat m `at'` (i * (xRow m) + j * (xCol m))
 {-# INLINE atM' #-}
 
 ------------------------------------------------------------------
 
+matrixFromVector :: Storable t => MatrixOrder -> Int -> Int -> Vector t -> Matrix t
 matrixFromVector _ 1 _ v@(dim->d) = Matrix { irows = 1, icols = d, xdat = v, xRow = d, xCol = 1 }
 matrixFromVector _ _ 1 v@(dim->d) = Matrix { irows = d, icols = 1, xdat = v, xRow = 1, xCol = d }
 matrixFromVector o r c v
@@ -388,18 +399,21 @@ subMatrix (r0,c0) (rt,ct) m
 
 --------------------------------------------------------------------------
 
+maxZ :: (Num t1, Ord t1, Foldable t) => t t1 -> t1
 maxZ xs = if minimum xs == 0 then 0 else maximum xs
 
+conformMs :: Element t => [Matrix t] -> [Matrix t]
 conformMs ms = map (conformMTo (r,c)) ms
   where
     r = maxZ (map rows ms)
     c = maxZ (map cols ms)
 
-
+conformVs :: Element t => [Vector t] -> [Vector t]
 conformVs vs = map (conformVTo n) vs
   where
     n = maxZ (map dim vs)
 
+conformMTo :: Element t => (Int, Int) -> Matrix t -> Matrix t
 conformMTo (r,c) m
     | size m == (r,c) = m
     | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
@@ -407,18 +421,24 @@ conformMTo (r,c) m
     | size m == (1,c) = repRows r m
     | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c)
 
+conformVTo :: Element t => Int -> Vector t -> Vector t
 conformVTo n v
     | dim v == n = v
     | dim v == 1 = constantD (v@>0) n
     | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
 
+repRows :: Element t => Int -> Matrix t -> Matrix t
 repRows n x = fromRows (replicate n (flatten x))
+repCols :: Element t => Int -> Matrix t -> Matrix t
 repCols n x = fromColumns (replicate n (flatten x))
 
+shSize :: Matrix t -> [Char]
 shSize = shDim . size
 
+shDim :: (Show a, Show a1) => (a1, a) -> [Char]
 shDim (r,c) = "(" ++ show r ++"x"++ show c ++")"
 
+emptyM :: Storable t => Int -> Int -> Matrix t
 emptyM r c = matrixFromVector RowMajor r c (fromList[])
 
 ----------------------------------------------------------------------
@@ -433,6 +453,11 @@ instance (Storable t, NFData t) => NFData (Matrix t)
 
 ---------------------------------------------------------------
 
+extractAux :: (Eq t3, Eq t2, TransArray c, Storable a, Storable t1,
+                Storable t, Num t3, Num t2, Integral t1, Integral t)
+           => (t3 -> t2 -> CInt -> Ptr t1 -> CInt -> Ptr t
+                  -> Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt))
+           -> MatrixOrder -> c -> t3 -> Vector t1 -> t2 -> Vector t -> IO (Matrix a)
 extractAux f ord m moder vr modec vc = do
     let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
         nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
@@ -452,6 +477,9 @@ foreign import ccall unsafe "extractL" c_extractL :: Extr Z
 
 ---------------------------------------------------------------
 
+setRectAux :: (TransArray c1, TransArray c)
+           => (CInt -> CInt -> Trans c1 (Trans c (IO CInt)))
+           -> Int -> Int -> c1 -> c -> IO ()
 setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"
 
 type SetRect x = I -> I -> x ::> x::> Ok
@@ -465,19 +493,29 @@ foreign import ccall unsafe "setRectL" c_setRectL :: SetRect Z
 
 --------------------------------------------------------------------------------
 
+sortG :: (Storable t, Storable a)
+      => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
 sortG f v = unsafePerformIO $ do
     r <- createVector (dim v)
     (v #! r) f #|"sortG"
     return r
 
+sortIdxD :: Vector Double -> Vector CInt
 sortIdxD = sortG c_sort_indexD
+sortIdxF :: Vector Float -> Vector CInt
 sortIdxF = sortG c_sort_indexF
+sortIdxI :: Vector CInt -> Vector CInt
 sortIdxI = sortG c_sort_indexI
+sortIdxL :: Vector Z -> Vector I
 sortIdxL = sortG c_sort_indexL
 
+sortValD :: Vector Double -> Vector Double
 sortValD = sortG c_sort_valD
+sortValF :: Vector Float -> Vector Float
 sortValF = sortG c_sort_valF
+sortValI :: Vector CInt -> Vector CInt
 sortValI = sortG c_sort_valI
+sortValL :: Vector Z -> Vector Z
 sortValL = sortG c_sort_valL
 
 foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt))
@@ -492,14 +530,21 @@ foreign import ccall unsafe "sort_valuesL" c_sort_valL :: Z :> Z :> Ok
 
 --------------------------------------------------------------------------------
 
+compareG :: (TransArray c, Storable t, Storable a)
+         => Trans c (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt)
+         -> c -> Vector t -> Vector a
 compareG f u v = unsafePerformIO $ do
     r <- createVector (dim v)
     (u # v #! r) f #|"compareG"
     return r
 
+compareD :: Vector Double -> Vector Double -> Vector CInt
 compareD = compareG c_compareD
+compareF :: Vector Float -> Vector Float -> Vector CInt
 compareF = compareG c_compareF
+compareI :: Vector CInt -> Vector CInt -> Vector CInt
 compareI = compareG c_compareI
+compareL :: Vector Z -> Vector Z -> Vector CInt
 compareL = compareG c_compareL
 
 foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt)))
@@ -509,16 +554,33 @@ foreign import ccall unsafe "compareL" c_compareL :: Z :> Z :> I :> Ok
 
 --------------------------------------------------------------------------------
 
+selectG :: (TransArray c, TransArray c1, TransArray c2, Storable t, Storable a)
+        => Trans c2 (Trans c1 (CInt -> Ptr t -> Trans c (CInt -> Ptr a -> IO CInt)))
+        -> c2 -> c1 -> Vector t -> c -> Vector a
 selectG f c u v w = unsafePerformIO $ do
     r <- createVector (dim v)
     (c # u # v # w #! r) f #|"selectG"
     return r
 
+selectD :: Vector CInt -> Vector Double -> Vector Double -> Vector Double -> Vector Double
 selectD = selectG c_selectD
+selectF :: Vector CInt -> Vector Float -> Vector Float -> Vector Float -> Vector Float
 selectF = selectG c_selectF
+selectI :: Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt
 selectI = selectG c_selectI
+selectL :: Vector CInt -> Vector Z -> Vector Z -> Vector Z -> Vector Z
 selectL = selectG c_selectL
+selectC :: Vector CInt
+        -> Vector (Complex Double)
+        -> Vector (Complex Double)
+        -> Vector (Complex Double)
+        -> Vector (Complex Double)
 selectC = selectG c_selectC
+selectQ :: Vector CInt
+        -> Vector (Complex Float)
+        -> Vector (Complex Float)
+        -> Vector (Complex Float)
+        -> Vector (Complex Float)
 selectQ = selectG c_selectQ
 
 type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt)))))
@@ -532,16 +594,29 @@ foreign import ccall unsafe "chooseL" c_selectL :: Sel Z
 
 ---------------------------------------------------------------------------
 
+remapG :: (TransArray c, TransArray c1, Storable t, Storable a)
+       => (CInt -> CInt -> CInt -> CInt -> Ptr t
+                -> Trans c1 (Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)))
+       -> Matrix t -> c1 -> c -> Matrix a
 remapG f i j m = unsafePerformIO $ do
     r <- createMatrix RowMajor (rows i) (cols i)
     (i # j # m #! r) f #|"remapG"
     return r
 
+remapD :: Matrix CInt -> Matrix CInt -> Matrix Double -> Matrix Double
 remapD = remapG c_remapD
+remapF :: Matrix CInt -> Matrix CInt -> Matrix Float -> Matrix Float
 remapF = remapG c_remapF
+remapI :: Matrix CInt -> Matrix CInt -> Matrix CInt -> Matrix CInt
 remapI = remapG c_remapI
+remapL :: Matrix CInt -> Matrix CInt -> Matrix Z -> Matrix Z
 remapL = remapG c_remapL
+remapC :: Matrix CInt
+       -> Matrix CInt
+       -> Matrix (Complex Double)
+       -> Matrix (Complex Double)
 remapC = remapG c_remapC
+remapQ :: Matrix CInt -> Matrix CInt -> Matrix (Complex Float) -> Matrix (Complex Float)
 remapQ = remapG c_remapQ
 
 type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt))))
@@ -555,6 +630,9 @@ foreign import ccall unsafe "remapL" c_remapL :: Rem Z
 
 --------------------------------------------------------------------------------
 
+rowOpAux :: (TransArray c, Storable a) =>
+            (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Trans c (IO CInt))
+         -> Int -> a -> Int -> Int -> Int -> Int -> c -> IO ()
 rowOpAux f c x i1 i2 j1 j2 m = do
     px <- newArray [x]
     (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp"
@@ -573,6 +651,9 @@ foreign import ccall unsafe "rowop_mod_int64_t" c_rowOpML :: Z -> RowOp Z
 
 --------------------------------------------------------------------------------
 
+gemmg :: (TransArray c1, TransArray c, TransArray c2, TransArray c3)
+      => Trans c3 (Trans c2 (Trans c1 (Trans c (IO CInt))))
+      -> c3 -> c2 -> c1 -> c -> IO ()
 gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"
 
 type Tgemm x = x :> x ::> x ::> x ::> Ok
@@ -588,6 +669,10 @@ foreign import ccall unsafe "gemm_mod_int64_t" c_gemmML :: Z -> Tgemm Z
 
 --------------------------------------------------------------------------------
 
+reorderAux :: (TransArray c, Storable t, Storable a1, Storable t1, Storable a) =>
+              (CInt -> Ptr a -> CInt -> Ptr t1
+                    -> Trans c (CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt))
+           -> Vector t1 -> c -> Vector t -> Vector a1
 reorderAux f s d v = unsafePerformIO $ do
     k <- createVector (dim s)
     r <- createVector (dim v)
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 9d51444..eb0c5a8 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -13,6 +13,9 @@
 {-# LANGUAGE TypeFamilies  #-}
 {-# LANGUAGE TypeOperators #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-missing-methods #-}
+
 {- |
 Module      :  Internal.Modular
 Copyright   :  (c) Alberto Ruiz 2015
diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs
index c9ef0c5..fd0a217 100644
--- a/packages/base/src/Internal/Numeric.hs
+++ b/packages/base/src/Internal/Numeric.hs
@@ -5,6 +5,8 @@
 {-# LANGUAGE FunctionalDependencies #-}
 {-# LANGUAGE UndecidableInstances #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Numeric
@@ -788,13 +790,7 @@ type instance RealOf (Complex Float) = Float
 type instance RealOf I = I
 type instance RealOf Z = Z
 
-type family ComplexOf x
-
-type instance ComplexOf Double = Complex Double
-type instance ComplexOf (Complex Double) = Complex Double
-
-type instance ComplexOf Float = Complex Float
-type instance ComplexOf (Complex Float) = Complex Float
+type ComplexOf x = Complex (RealOf x)
 
 type family SingleOf x
 
diff --git a/packages/base/src/Internal/ST.hs b/packages/base/src/Internal/ST.hs
index 544c9e4..7d54e6d 100644
--- a/packages/base/src/Internal/ST.hs
+++ b/packages/base/src/Internal/ST.hs
@@ -81,6 +81,8 @@ unsafeFreezeVector :: (Storable t) => STVector s t -> ST s (Vector t)
 unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
 
 {-# INLINE safeIndexV #-}
+safeIndexV :: Storable t2
+           => (STVector s t2 -> Int -> t) -> STVector t1 t2 -> Int -> t
 safeIndexV f (STVector v) k
     | k < 0 || k>= dim v = error $ "out of range error in vector (dim="
                                    ++show (dim v)++", pos="++show k++")"
@@ -150,9 +152,12 @@ unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
 freezeMatrix :: (Element t) => STMatrix s t -> ST s (Matrix t)
 freezeMatrix m = liftSTMatrix id m
 
+cloneMatrix :: Element t => Matrix t -> IO (Matrix t)
 cloneMatrix m = copy (orderOf m) m
 
 {-# INLINE safeIndexM #-}
+safeIndexM :: (STMatrix s t2 -> Int -> Int -> t)
+           -> STMatrix t1 t2 -> Int -> Int -> t
 safeIndexM f (STMatrix m) r c
     | r<0 || r>=rows m ||
       c<0 || c>=cols m = error $ "out of range error in matrix (size="
@@ -184,6 +189,7 @@ data ColRange = AllCols
               | Col Int
               | FromCol Int
 
+getColRange :: Int -> ColRange -> (Int, Int)
 getColRange c AllCols = (0,c-1)
 getColRange c (ColRange a b) = (a `mod` c, b `mod` c)
 getColRange c (Col a) = (a `mod` c, a `mod` c)
@@ -194,6 +200,7 @@ data RowRange = AllRows
               | Row Int
               | FromRow Int
 
+getRowRange :: Int -> RowRange -> (Int, Int)
 getRowRange r AllRows = (0,r-1)
 getRowRange r (RowRange a b) = (a `mod` r, b `mod` r)
 getRowRange r (Row a) = (a `mod` r, a `mod` r)
@@ -223,6 +230,7 @@ rowOper (SWAP i1 i2 r) (STMatrix m) = unsafeIOToST $ rowOp 2 0 i1' i2' j1 j2 m
     i2' = i2 `mod` (rows m)
 
 
+extractMatrix :: Element a => STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
 extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[i1,i2]) 0 (idxs[j1,j2]))
   where
     (i1,i2) = getRowRange (rows m) rr
@@ -231,6 +239,7 @@ extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[
 -- | r0 c0 height width
 data Slice s t = Slice (STMatrix s t) Int Int Int Int
 
+slice :: Element a => Slice t a -> Matrix a
 slice (Slice (STMatrix m) r0 c0 nr nc) = subMatrix (r0,c0) (nr,nc) m
 
 gemmm :: Element t => t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
@@ -238,7 +247,7 @@ gemmm beta (slice->r) alpha (slice->a) (slice->b) = res
   where
     res = unsafeIOToST (gemm v a b r)
     v = fromList [alpha,beta]
-    
+
 
 mutable :: Element t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u)
 mutable f a = runST $ do
@@ -246,4 +255,3 @@ mutable f a = runST $ do
    info <- f (rows a, cols a) x
    r <- unsafeFreezeMatrix x
    return (r,info)
-
diff --git a/packages/base/src/Internal/Sparse.hs b/packages/base/src/Internal/Sparse.hs
index a8a5fe0..fbea11a 100644
--- a/packages/base/src/Internal/Sparse.hs
+++ b/packages/base/src/Internal/Sparse.hs
@@ -2,6 +2,8 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FlexibleInstances #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 module Internal.Sparse(
     GMatrix(..), CSR(..), mkCSR, fromCSR,
     mkSparse, mkDiagR, mkDense,
diff --git a/packages/base/src/Internal/Static.hs b/packages/base/src/Internal/Static.hs
index 357645e..566506c 100644
--- a/packages/base/src/Internal/Static.hs
+++ b/packages/base/src/Internal/Static.hs
@@ -15,6 +15,9 @@
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE DeriveGeneric #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}
+
 {- |
 Module      :  Internal.Static
 Copyright   :  (c) Alberto Ruiz 2006-14
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 959e58f..f642e8d 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -6,6 +6,8 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE ViewPatterns #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
 
 -----------------------------------------------------------------------------
 {- |
diff --git a/packages/base/src/Internal/Vector.hs b/packages/base/src/Internal/Vector.hs
index e1e4aa8..6271bb6 100644
--- a/packages/base/src/Internal/Vector.hs
+++ b/packages/base/src/Internal/Vector.hs
@@ -1,6 +1,7 @@
 {-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-}
 {-# LANGUAGE TypeSynonymInstances #-}
 
+{-# OPTIONS_GHC -fno-warn-orphans #-}
 
 -- |
 -- Module      :  Internal.Vector
@@ -40,6 +41,7 @@ import qualified Data.Vector.Storable as Vector
 import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith)
 
 import Data.Binary
+import Data.Binary.Put
 import Control.Monad(replicateM)
 import qualified Data.ByteString.Internal as BS
 import Data.Vector.Storable.Internal(updPtr)
@@ -92,6 +94,7 @@ createVector n = do
 
 -}
 
+safeRead :: Storable a => Vector a -> (Ptr a -> IO c) -> c
 safeRead v = inlinePerformIO . unsafeWith v
 {-# INLINE safeRead #-}
 
@@ -287,11 +290,13 @@ foldVectorWithIndex f x v = unsafePerformIO $
         go (dim v -1) x
 {-# INLINE foldVectorWithIndex #-}
 
+foldLoop :: (Int -> t -> t) -> t -> Int -> t
 foldLoop f s0 d = go (d - 1) s0
      where
        go 0 s = f (0::Int) s
        go !j !s = go (j - 1) (f j s)
 
+foldVectorG :: Storable t1 => (Int -> (Int -> t1) -> t -> t) -> t -> Vector t1 -> t
 foldVectorG f s0 v = foldLoop g s0 (dim v)
     where g !k !s = f k (safeRead v . flip peekElemOff) s
           {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
@@ -394,8 +399,10 @@ chunks d = let c = d `div` chunk
                m = d `mod` chunk
            in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk)
 
+putVector :: (Storable t, Binary t) => Vector t -> Data.Binary.Put.PutM ()
 putVector v = mapM_ put $! toList v
 
+getVector :: (Storable a, Binary a) => Int -> Get (Vector a)
 getVector d = do
               xs <- replicateM d get
               return $! fromList xs
diff --git a/packages/base/src/Internal/Vectorized.hs b/packages/base/src/Internal/Vectorized.hs
index 2990173..32430c6 100644
--- a/packages/base/src/Internal/Vectorized.hs
+++ b/packages/base/src/Internal/Vectorized.hs
@@ -28,12 +28,15 @@ import System.IO.Unsafe(unsafePerformIO)
 import Control.Monad(when)
 
 infixr 1 #
+(#) :: TransArray c => c -> (b -> IO r) -> TransRaw c b -> IO r
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
+(#!) :: (TransArray c, TransArray c1) => c1 -> c -> TransRaw c1 (TransRaw c (IO r)) -> IO r
 a #! b = a # b # id
 {-# INLINE (#!) #-}
 
+fromei :: Enum a => a -> CInt
 fromei x = fromIntegral (fromEnum x) :: CInt
 
 data FunCodeV = Sin
@@ -100,10 +103,20 @@ sumQ = sumg c_sumQ
 sumC :: Vector (Complex Double) -> Complex Double
 sumC = sumg c_sumC
 
+sumI :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr I -> I :> Ok)
+        , TransArray c
+        , Storable a
+        )
+     => I -> c -> a
 sumI m = sumg (c_sumI m)
 
+sumL :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr Z -> Z :> Ok)
+        , TransArray c
+        , Storable a
+        ) => Z -> c -> a
 sumL m = sumg (c_sumL m)
 
+sumg :: (TransArray c, Storable a) => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a
 sumg f x = unsafePerformIO $ do
     r <- createVector 1
     (x #! r) f #| "sum"
@@ -140,6 +153,8 @@ prodI = prodg . c_prodI
 prodL :: Z-> Vector Z -> Z
 prodL = prodg . c_prodL
 
+prodg :: (TransArray c, Storable a)
+      => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a
 prodg f x = unsafePerformIO $ do
     r <- createVector 1
     (x #! r) f #| "prod"
@@ -155,16 +170,25 @@ foreign import ccall unsafe "prodL" c_prodL :: Z -> TVV Z
 
 ------------------------------------------------------------------
 
+toScalarAux :: (Enum a, TransArray c, Storable a1)
+            => (CInt -> TransRaw c (CInt -> Ptr a1 -> IO CInt)) -> a -> c -> a1
 toScalarAux fun code v = unsafePerformIO $ do
     r <- createVector 1
     (v #! r) (fun (fromei code)) #|"toScalarAux"
     return (r @> 0)
 
+
+vectorMapAux :: (Enum a, Storable t, Storable a1)
+             => (CInt -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)
+             -> a -> Vector t -> Vector a1
 vectorMapAux fun code v = unsafePerformIO $ do
     r <- createVector (dim v)
     (v #! r) (fun (fromei code)) #|"vectorMapAux"
     return r
 
+vectorMapValAux :: (Enum a, Storable a2, Storable t, Storable a1)
+                => (CInt -> Ptr a2 -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)
+                -> a -> a2 -> Vector t -> Vector a1
 vectorMapValAux fun code val v = unsafePerformIO $ do
     r <- createVector (dim v)
     pval <- newArray [val]
@@ -172,6 +196,9 @@ vectorMapValAux fun code val v = unsafePerformIO $ do
     free pval
     return r
 
+vectorZipAux :: (Enum a, TransArray c, Storable t, Storable a1)
+             => (CInt -> CInt -> Ptr t -> TransRaw c (CInt -> Ptr a1 -> IO CInt))
+             -> a -> Vector t -> c -> Vector a1
 vectorZipAux fun code u v = unsafePerformIO $ do
     r <- createVector (dim u)
     (u # v #! r) (fun (fromei code)) #|"vectorZipAux"
@@ -378,6 +405,7 @@ foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> D
 
 --------------------------------------------------------------------------------
 
+roundVector :: Vector Double -> Vector Double
 roundVector v = unsafePerformIO $ do
     r <- createVector (dim v)
     (v #! r) c_round_vector #|"roundVector"
@@ -433,6 +461,8 @@ long2intV :: Vector Z -> Vector I
 long2intV = tog c_long2int
 
 
+tog :: (Storable t, Storable a)
+    => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
 tog f v = unsafePerformIO $ do
     r <- createVector (dim v)
     (v #! r) f #|"tog"
@@ -452,6 +482,8 @@ foreign import ccall unsafe "long2int"    c_long2int    :: Z :> I :> Ok
 
 ---------------------------------------------------------------
 
+stepg :: (Storable t, Storable a)
+      => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
 stepg f v = unsafePerformIO $ do
     r <- createVector (dim v)
     (v #! r) f #|"step"
@@ -477,6 +509,8 @@ foreign import ccall unsafe "stepL" c_stepL :: TVV Z
 
 --------------------------------------------------------------------------------
 
+conjugateAux :: (Storable t, Storable a)
+             => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
 conjugateAux fun x = unsafePerformIO $ do
     v <- createVector (dim x)
     (x #! v) fun #|"conjugateAux"
@@ -502,6 +536,8 @@ cloneVector v = do
 
 --------------------------------------------------------------------------------
 
+constantAux :: (Storable a1, Storable a)
+            => (Ptr a1 -> CInt -> Ptr a -> IO CInt) -> a1 -> Int -> Vector a
 constantAux fun x n = unsafePerformIO $ do
     v <- createVector n
     px <- newArray [x]
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 520eeb7..91923e9 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -1,6 +1,8 @@
 {-# LANGUAGE CPP #-}
 {-# LANGUAGE FlexibleContexts #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
index 3a84645..57e5cf1 100644
--- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -28,7 +28,9 @@ infixr 8 <·>
 (<·>) :: Numeric t => Vector t -> Vector t -> t
 (<·>) = dot
 
+app :: Numeric t => Matrix t -> Vector t -> Vector t
 app m v = m #> v
 
+mul :: Numeric t => Matrix t -> Matrix t -> Matrix t
 mul a b = a <> b
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index e328904..2e05c90 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -14,6 +14,8 @@
 {-# LANGUAGE GADTs #-}
 {-# LANGUAGE TypeFamilies #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
 
 {- |
 Module      :  Numeric.LinearAlgebra.Static
diff --git a/packages/base/src/Numeric/Matrix.hs b/packages/base/src/Numeric/Matrix.hs
index 06da150..6e3db61 100644
--- a/packages/base/src/Numeric/Matrix.hs
+++ b/packages/base/src/Numeric/Matrix.hs
@@ -4,6 +4,8 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Matrix
@@ -35,6 +37,7 @@ import Data.List(partition)
 import qualified Data.Foldable as F
 import qualified Data.Semigroup as S
 import Internal.Chain
+import Foreign.Storable(Storable)
 
 
 -------------------------------------------------------------------
@@ -80,8 +83,16 @@ instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matri
 
 --------------------------------------------------------------------------------
 
+isScalar :: Matrix t -> Bool
 isScalar m = rows m == 1 && cols m == 1
 
+adaptScalarM :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2)
+             => (t1 -> Matrix t2 -> t)
+             -> (Matrix t1 -> Matrix t2 -> t)
+             -> (Matrix t1 -> t2 -> t)
+             -> Matrix t1
+             -> Matrix t2
+             -> t
 adaptScalarM f1 f2 f3 x y
     | isScalar x = f1   (x @@>(0,0) ) y
     | isScalar y = f3 x (y @@>(0,0) )
@@ -96,7 +107,7 @@ instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matr
   where
     mempty = 1
     mappend = adaptScalarM scale mXm (flip scale)
-    
+
     mconcat xs = work (partition isScalar xs)
       where
         work (ss,[]) = product ss
@@ -106,4 +117,3 @@ instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matr
             | otherwise              = scale x00 m
           where
             x00 = x @@> (0,0)
-
diff --git a/packages/base/src/Numeric/Vector.hs b/packages/base/src/Numeric/Vector.hs
index 017196c..1e5877d 100644
--- a/packages/base/src/Numeric/Vector.hs
+++ b/packages/base/src/Numeric/Vector.hs
@@ -3,6 +3,9 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Vector
@@ -14,7 +17,7 @@
 --
 -- Provides instances of standard classes 'Show', 'Read', 'Eq',
 -- 'Num', 'Fractional',  and 'Floating' for 'Vector'.
--- 
+--
 -----------------------------------------------------------------------------
 
 module Numeric.Vector () where
@@ -23,9 +26,17 @@ import Internal.Vectorized
 import Internal.Vector
 import Internal.Numeric
 import Internal.Conversion
+import Foreign.Storable(Storable)
 
 -------------------------------------------------------------------
 
+adaptScalar :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2)
+            => (t1 -> Vector t2 -> t)
+            -> (Vector t1 -> Vector t2 -> t)
+            -> (Vector t1 -> t2 -> t)
+            -> Vector t1
+            -> Vector t2
+            -> t
 adaptScalar f1 f2 f3 x y
     | dim x == 1 = f1   (x@>0) y
     | dim y == 1 = f3 x (y@>0)
@@ -172,4 +183,3 @@ instance Floating (Vector (Complex Float)) where
     sqrt  = vectorMapQ Sqrt
     (**)  = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
     pi    = fromList [pi]
-
diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal
index 6b0032b..ca93775 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -1,6 +1,6 @@
 Name:               hmatrix-glpk
-Version:            0.6.0.0
-License:            GPL
+Version:            0.19.0.0
+License:            GPL-3
 License-file:       LICENSE
 Author:             Alberto Ruiz
 Maintainer:         Alberto Ruiz <aruiz@um.es>
diff --git a/packages/gsl/CHANGELOG b/packages/gsl/CHANGELOG
index 091dc0e..a2fe038 100644
--- a/packages/gsl/CHANGELOG
+++ b/packages/gsl/CHANGELOG
@@ -2,7 +2,7 @@
 --------
 
     * Added interpolation modules
-    
+
     * Added simulated annealing module
 
     * Added odeSolveVWith
@@ -11,4 +11,3 @@
 --------
 
     * The modules Numeric.GSL.* have been moved from hmatrix to the new package hmatrix-gsl.
-
diff --git a/packages/gsl/hmatrix-gsl.cabal b/packages/gsl/hmatrix-gsl.cabal
index d463ee8..76db835 100644
--- a/packages/gsl/hmatrix-gsl.cabal
+++ b/packages/gsl/hmatrix-gsl.cabal
@@ -1,21 +1,18 @@
 Name:               hmatrix-gsl
-Version:            0.18.0.1
-License:            GPL
-License-file:       LICENSE
-Author:             Alberto Ruiz
-Maintainer:         Alberto Ruiz <aruiz@um.es>
-Stability:          provisional
-Homepage:           https://github.com/albertoruiz/hmatrix
+Version:            0.19.0.0
 Synopsis:           Numerical computation
 Description:        Purely functional interface to selected numerical computations,
                     internally implemented using GSL.
-
+Homepage:           https://github.com/albertoruiz/hmatrix
+license:            GPL-3
+license-file:       LICENSE
+Author:             Alberto Ruiz
+Maintainer:         Alberto Ruiz <aruiz@um.es>
+Stability:          provisional
 Category:           Math
-tested-with:        GHC ==7.8
-
-cabal-version:      >=1.8
-
 build-type:         Simple
+cabal-version:      >=1.18
+
 
 extra-source-files: src/Numeric/GSL/gsl-ode.c
 
@@ -33,9 +30,6 @@ library
     Build-Depends:      base<5, hmatrix>=0.18, array, vector,
                         process, random
 
-
-    Extensions:         ForeignFunctionInterface
-
     hs-source-dirs:     src
     Exposed-modules:    Numeric.GSL.Differentiation,
                         Numeric.GSL.Integration,
@@ -98,6 +92,8 @@ library
     else
       pkgconfig-depends: gsl
 
+    default-language:    Haskell2010
+
 
 source-repository head
     type:     git
diff --git a/packages/gsl/src/Graphics/Plot.hs b/packages/gsl/src/Graphics/Plot.hs
index d2ea192..e422912 100644
--- a/packages/gsl/src/Graphics/Plot.hs
+++ b/packages/gsl/src/Graphics/Plot.hs
@@ -1,3 +1,5 @@
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Graphics.Plot
diff --git a/packages/gsl/src/Numeric/GSL/Fitting.hs b/packages/gsl/src/Numeric/GSL/Fitting.hs
index 8f2eae3..a732c25 100644
--- a/packages/gsl/src/Numeric/GSL/Fitting.hs
+++ b/packages/gsl/src/Numeric/GSL/Fitting.hs
@@ -1,6 +1,8 @@
 {-# LANGUAGE CPP #-}
 {-# LANGUAGE FlexibleContexts #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 {- |
 Module      :  Numeric.GSL.Fitting
 Copyright   :  (c) Alberto Ruiz 2010
diff --git a/packages/gsl/src/Numeric/GSL/Fourier.hs b/packages/gsl/src/Numeric/GSL/Fourier.hs
index bffab87..ed7353a 100644
--- a/packages/gsl/src/Numeric/GSL/Fourier.hs
+++ b/packages/gsl/src/Numeric/GSL/Fourier.hs
@@ -1,5 +1,7 @@
 {-# LANGUAGE TypeFamilies #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 {- |
 Module      : Numeric.GSL.Fourier
 Copyright   :  (c) Alberto Ruiz 2006
diff --git a/packages/gsl/src/Numeric/GSL/Integration.hs b/packages/gsl/src/Numeric/GSL/Integration.hs
index 9c1d43a..0a1b4c6 100644
--- a/packages/gsl/src/Numeric/GSL/Integration.hs
+++ b/packages/gsl/src/Numeric/GSL/Integration.hs
@@ -1,3 +1,5 @@
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 {- |
 Module      :  Numeric.GSL.Integration
 Copyright   :  (c) Alberto Ruiz 2006
diff --git a/packages/gsl/src/Numeric/GSL/Internal.hs b/packages/gsl/src/Numeric/GSL/Internal.hs
index f70e167..e1f8d95 100644
--- a/packages/gsl/src/Numeric/GSL/Internal.hs
+++ b/packages/gsl/src/Numeric/GSL/Internal.hs
@@ -1,5 +1,8 @@
 {-# LANGUAGE FlexibleContexts #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 -- |
 -- Module      :  Numeric.GSL.Internal
 -- Copyright   :  (c) Alberto Ruiz 2009
@@ -128,8 +131,7 @@ type TVM = TV (TM Res)
 ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
 
 vec x f = unsafeWith x $ \p -> do
-    let v g = do
-        g (fi $ V.length x) p
+    let v g = g (fi $ V.length x) p
     f v
 {-# INLINE vec #-}
 
diff --git a/packages/gsl/src/Numeric/GSL/Interpolation.hs b/packages/gsl/src/Numeric/GSL/Interpolation.hs
index 6f02405..484d2a2 100644
--- a/packages/gsl/src/Numeric/GSL/Interpolation.hs
+++ b/packages/gsl/src/Numeric/GSL/Interpolation.hs
@@ -1,5 +1,7 @@
 {-# LANGUAGE MagicHash, UnboxedTuples #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 {- |
 Module      :  Numeric.GSL.Interpolation
 Copyright   :  (c) Matthew Peddie 2015
diff --git a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
index 1bf357b..aee64f7 100644
--- a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
@@ -1,3 +1,6 @@
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.LinearAlgebra
diff --git a/packages/gsl/src/Numeric/GSL/Minimization.hs b/packages/gsl/src/Numeric/GSL/Minimization.hs
index a0e5306..1fd951b 100644
--- a/packages/gsl/src/Numeric/GSL/Minimization.hs
+++ b/packages/gsl/src/Numeric/GSL/Minimization.hs
@@ -1,3 +1,5 @@
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 {-# LANGUAGE FlexibleContexts #-}
 
 
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
index 987d47e..a1ccd38 100644
--- a/packages/gsl/src/Numeric/GSL/ODE.hs
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -1,5 +1,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
 
 {- |
 Module      :  Numeric.GSL.ODE
diff --git a/packages/gsl/src/Numeric/GSL/Root.hs b/packages/gsl/src/Numeric/GSL/Root.hs
index 724f32f..9cdb061 100644
--- a/packages/gsl/src/Numeric/GSL/Root.hs
+++ b/packages/gsl/src/Numeric/GSL/Root.hs
@@ -1,5 +1,7 @@
 {-# LANGUAGE  FlexibleContexts #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 {- |
 Module      :  Numeric.GSL.Root
 Copyright   :  (c) Alberto Ruiz 2009
diff --git a/packages/gsl/src/Numeric/GSL/Vector.hs b/packages/gsl/src/Numeric/GSL/Vector.hs
index b1c0106..2ca7cc0 100644
--- a/packages/gsl/src/Numeric/GSL/Vector.hs
+++ b/packages/gsl/src/Numeric/GSL/Vector.hs
@@ -1,3 +1,6 @@
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Vector
diff --git a/packages/sparse/hmatrix-sparse.cabal b/packages/sparse/hmatrix-sparse.cabal
index 55eb424..4399b72 100644
--- a/packages/sparse/hmatrix-sparse.cabal
+++ b/packages/sparse/hmatrix-sparse.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-sparse
-Version:            0.1.0
+Version:            0.19.0.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -19,7 +19,7 @@ build-type:         Simple
 
 
 library
-    Build-Depends:      base, hmatrix>=0.16
+    Build-Depends:      base<5, hmatrix>=0.16
 
     hs-source-dirs:     src
 
diff --git a/packages/special/hmatrix-special.cabal b/packages/special/hmatrix-special.cabal
index 2848e39..0890bc7 100644
--- a/packages/special/hmatrix-special.cabal
+++ b/packages/special/hmatrix-special.cabal
@@ -1,6 +1,6 @@
 Name:               hmatrix-special
-Version:            0.4.0.1
-License:            GPL
+Version:            0.19.0.0
+License:            GPL-3
 License-file:       LICENSE
 Author:             Alberto Ruiz
 Maintainer:         Alberto Ruiz <aruiz@um.es>
diff --git a/packages/special/lib/Numeric/GSL/Special/Bessel.hs b/packages/special/lib/Numeric/GSL/Special/Bessel.hs
index 70066f8..84d4cf5 100644
--- a/packages/special/lib/Numeric/GSL/Special/Bessel.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Bessel.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Bessel
diff --git a/packages/special/lib/Numeric/GSL/Special/Coulomb.hs b/packages/special/lib/Numeric/GSL/Special/Coulomb.hs
index 6904739..3bd3ed6 100644
--- a/packages/special/lib/Numeric/GSL/Special/Coulomb.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Coulomb.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Coulomb
diff --git a/packages/special/lib/Numeric/GSL/Special/Coupling.hs b/packages/special/lib/Numeric/GSL/Special/Coupling.hs
index ad120cc..e8d9aef 100644
--- a/packages/special/lib/Numeric/GSL/Special/Coupling.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Coupling.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Coupling
diff --git a/packages/special/lib/Numeric/GSL/Special/Exp.hs b/packages/special/lib/Numeric/GSL/Special/Exp.hs
index b6dfeef..54033c5 100644
--- a/packages/special/lib/Numeric/GSL/Special/Exp.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Exp.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Exp
diff --git a/packages/special/lib/Numeric/GSL/Special/Gamma.hs b/packages/special/lib/Numeric/GSL/Special/Gamma.hs
index 41e24f0..55950cc 100644
--- a/packages/special/lib/Numeric/GSL/Special/Gamma.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Gamma.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Gamma
diff --git a/packages/special/lib/Numeric/GSL/Special/Gegenbauer.hs b/packages/special/lib/Numeric/GSL/Special/Gegenbauer.hs
index fb8bf3f..1dae1f1 100644
--- a/packages/special/lib/Numeric/GSL/Special/Gegenbauer.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Gegenbauer.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Gegenbauer
diff --git a/packages/special/lib/Numeric/GSL/Special/Legendre.hs b/packages/special/lib/Numeric/GSL/Special/Legendre.hs
index 927fa2c..5f7d2b0 100644
--- a/packages/special/lib/Numeric/GSL/Special/Legendre.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Legendre.hs
@@ -1,4 +1,7 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Legendre
diff --git a/packages/special/lib/Numeric/GSL/Special/Trig.hs b/packages/special/lib/Numeric/GSL/Special/Trig.hs
index f2c1519..754bed1 100644
--- a/packages/special/lib/Numeric/GSL/Special/Trig.hs
+++ b/packages/special/lib/Numeric/GSL/Special/Trig.hs
@@ -1,4 +1,8 @@
 {-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
+
 ------------------------------------------------------------
 -- |
 -- Module      :  Numeric.GSL.Special.Trig
diff --git a/packages/sundials/ChangeLog.md b/packages/sundials/ChangeLog.md
new file mode 100644
index 0000000..7b15777
--- /dev/null
+++ b/packages/sundials/ChangeLog.md
@@ -0,0 +1,5 @@
+# Revision history for hmatrix-sundials
+
+## 0.1.0.0  -- 2018-04-21
+
+* First version. Released on an unsuspecting world. Just Runge-Kutta methods to start with.
diff --git a/packages/sundials/LICENSE b/packages/sundials/LICENSE
new file mode 100644
index 0000000..a162e98
--- /dev/null
+++ b/packages/sundials/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2018, Dominic Steinitz, Novadiscovery
+
+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 Dominic Steinitz 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/packages/sundials/README.md b/packages/sundials/README.md
new file mode 100644
index 0000000..2fac5c2
--- /dev/null
+++ b/packages/sundials/README.md
@@ -0,0 +1,8 @@
+Currently only an interface to the Runge-Kutta methods:
+[ARKode](https://computation.llnl.gov/projects/sundials/arkode)
+
+The interface is almost certainly going to change. Sundials gives a
+rich set of "combinators" for controlling the solution of your problem
+and reporting on how it performed. The idea is to initially mimic
+hmatrix-gsl and add extra, richer functions but ultimately upgrade the
+whole interface both for sundials and for gsl.
diff --git a/packages/sundials/Setup.hs b/packages/sundials/Setup.hs
new file mode 100644
index 0000000..9a994af
--- /dev/null
+++ b/packages/sundials/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/packages/sundials/diagrams/brusselator.png b/packages/sundials/diagrams/brusselator.png
new file mode 100644
index 0000000..740cacb
--- /dev/null
+++ b/packages/sundials/diagrams/brusselator.png
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
new file mode 100644
index 0000000..cd2be4e
--- /dev/null
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -0,0 +1,61 @@
+name:                hmatrix-sundials
+version:             0.19.0.0
+synopsis:            hmatrix interface to sundials
+description:         An interface to the solving suite SUNDIALS. Currently, it
+                     mimics the solving interface in hmstrix-gsl but
+                     provides more diagnostic information and the
+                     Butcher Tableaux (for Runge-Kutta methods).
+homepage:            https://github.com/idontgetoutmuch/hmatrix/tree/sundials
+license:             BSD3
+license-file:        LICENSE
+author:              Dominic Steinitz
+maintainer:          dominic@steinitz.org
+copyright:           Dominic Steinitz 2018, Novadiscovery 2018
+category:            Math
+build-type:          Simple
+extra-source-files:  ChangeLog.md, README.md, diagrams/*.png
+extra-doc-files:     diagrams/*.png
+cabal-version:       >=1.18
+
+
+library
+  build-depends:       base >=4.10 && <4.11,
+                       inline-c >=0.6 && <0.7,
+                       vector >=0.12 && <0.13,
+                       template-haskell >=2.12 && <2.13,
+                       containers >=0.5 && <0.6,
+                       hmatrix>=0.18
+  extra-libraries:     sundials_arkode,
+                       sundials_cvode
+  other-extensions:    QuasiQuotes
+  hs-source-dirs:      src
+  exposed-modules:     Numeric.Sundials.ODEOpts,
+                       Numeric.Sundials.ARKode.ODE,
+                       Numeric.Sundials.CVode.ODE
+  other-modules:       Numeric.Sundials.Arkode
+  c-sources:           src/helpers.c src/helpers.h
+  default-language:    Haskell2010
+
+test-suite hmatrix-sundials-testsuite
+  type:                exitcode-stdio-1.0
+  main-is:             Main.hs
+  other-modules:       Numeric.Sundials.ODEOpts,
+                       Numeric.Sundials.ARKode.ODE,
+                       Numeric.Sundials.CVode.ODE,
+                       Numeric.Sundials.Arkode
+  build-depends:       base >=4.10 && <4.11,
+                       inline-c >=0.6 && <0.7,
+                       vector >=0.12 && <0.13,
+                       template-haskell >=2.12 && <2.13,
+                       containers >=0.5 && <0.6,
+                       hmatrix>=0.18,
+                       plots,
+                       diagrams-lib,
+                       diagrams-rasterific,
+                       lens,
+                       hspec
+  hs-source-dirs:      src
+  extra-libraries:     sundials_arkode,
+                       sundials_cvode
+  c-sources:           src/helpers.c src/helpers.h
+  default-language:    Haskell2010
diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
new file mode 100644
index 0000000..16c21c5
--- /dev/null
+++ b/packages/sundials/src/Main.hs
@@ -0,0 +1,186 @@
+{-# OPTIONS_GHC -Wall #-}
+
+import qualified Numeric.Sundials.ARKode.ODE as ARK
+import qualified Numeric.Sundials.CVode.ODE  as CV
+import           Numeric.LinearAlgebra
+
+import           Plots as P
+import qualified Diagrams.Prelude as D
+import           Diagrams.Backend.Rasterific
+
+import           Control.Lens
+
+import           Test.Hspec
+
+
+lorenz :: Double -> [Double] -> [Double]
+lorenz _t u = [ sigma * (y - x)
+              , x * (rho - z) - y
+              , x * y - beta * z
+              ]
+  where
+    rho = 28.0
+    sigma = 10.0
+    beta = 8.0 / 3.0
+    x = u !! 0
+    y = u !! 1
+    z = u !! 2
+
+_lorenzJac :: Double -> Vector Double -> Matrix Double
+_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y
+                        , sigma   , -1.0   , x
+                        , 0.0     , (-x)   , (-beta)
+                        ]
+  where
+    rho = 28.0
+    sigma = 10.0
+    beta = 8.0 / 3.0
+    x = u ! 0
+    y = u ! 1
+    z = u ! 2
+
+brusselator :: Double -> [Double] -> [Double]
+brusselator _t x = [ a - (w + 1) * u + v * u * u
+                   , w * u - v * u * u
+                   , (b - w) / eps - w * u
+                   ]
+  where
+    a = 1.0
+    b = 3.5
+    eps = 5.0e-6
+    u = x !! 0
+    v = x !! 1
+    w = x !! 2
+
+_brussJac :: Double -> Vector Double -> Matrix Double
+_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
+                       , u * u                     , (-(u * u))     , 0.0
+                       , (-u)                      , u              , (-1.0) / eps - u
+                       ]
+  where
+    y = toList x
+    u = y !! 0
+    v = y !! 1
+    w = y !! 2
+    eps = 5.0e-6
+
+stiffish :: Double -> [Double] -> [Double]
+stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
+  where
+    lamda = -100.0
+    u = v !! 0
+
+stiffishV :: Double -> Vector Double -> Vector Double
+stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
+  where
+    lamda = -100.0
+    u = v ! 0
+
+_stiffJac :: Double -> Vector Double -> Matrix Double
+_stiffJac _t _v = (1><1) [ lamda ]
+  where
+    lamda = -100.0
+
+predatorPrey :: Double -> [Double] -> [Double]
+predatorPrey _t v = [ x * a - b * x * y
+                    , d * x * y - c * y - e * y * z
+                    , (-f) * z + g * y * z
+                    ]
+  where
+    x = v!!0
+    y = v!!1
+    z = v!!2
+    a = 1.0
+    b = 1.0
+    c = 1.0
+    d = 1.0
+    e = 1.0
+    f = 1.0
+    g = 1.0
+
+lSaxis :: [[Double]] -> P.Axis B D.V2 Double
+lSaxis xs = P.r2Axis &~ do
+  let ts = xs!!0
+      us = xs!!1
+      vs = xs!!2
+      ws = xs!!3
+  P.linePlot' $ zip ts us
+  P.linePlot' $ zip ts vs
+  P.linePlot' $ zip ts ws
+
+kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
+kSaxis xs = P.r2Axis &~ do
+  P.linePlot' xs
+
+main :: IO ()
+main = do
+
+  let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
+  renderRasterific "diagrams/brusselator.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
+
+  let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
+  renderRasterific "diagrams/brusselatorA.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))
+
+  let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
+  renderRasterific "diagrams/stiffish.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2))
+
+  let res2a = ARK.odeSolveV (ARK.SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
+
+  let res2b = ARK.odeSolveV (ARK.TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
+
+  let maxDiffA = maximum $ map abs $
+                 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)
+
+  let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
+
+  let maxDiffB = maximum $ map abs $
+                 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0)
+
+  let maxDiffC = maximum $ map abs $
+                 zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0)
+
+  let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
+
+  renderRasterific "diagrams/lorenz.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))
+
+  renderRasterific "diagrams/lorenz1.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2))
+
+  renderRasterific "diagrams/lorenz2.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2))
+
+  let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
+
+  renderRasterific "diagrams/predatorPrey.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1))
+
+  renderRasterific "diagrams/predatorPrey1.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2))
+
+  renderRasterific "diagrams/predatorPrey2.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2))
+
+  let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])
+
+  let maxDiffPpA = maximum $ map abs $
+                   zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0)
+
+  hspec $ describe "Compare results" $ do
+    it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6
+    it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6
+    it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6
+    it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3
+
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
new file mode 100644
index 0000000..fafc237
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -0,0 +1,903 @@
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE MultiWayIf #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Sundials.ARKode.ODE
+-- Copyright   :  Dominic Steinitz 2018,
+--                Novadiscovery 2018
+-- License     :  BSD
+-- Maintainer  :  Dominic Steinitz
+-- Stability   :  provisional
+--
+-- Solution of ordinary differential equation (ODE) initial value problems.
+-- See <https://computation.llnl.gov/projects/sundials/sundials-software> for more detail.
+--
+-- A simple example:
+--
+-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
+--
+-- @
+-- import           Numeric.Sundials.ARKode.ODE
+-- import           Numeric.LinearAlgebra
+--
+-- import           Plots as P
+-- import qualified Diagrams.Prelude as D
+-- import           Diagrams.Backend.Rasterific
+--
+-- brusselator :: Double -> [Double] -> [Double]
+-- brusselator _t x = [ a - (w + 1) * u + v * u * u
+--                    , w * u - v * u * u
+--                    , (b - w) / eps - w * u
+--                    ]
+--   where
+--     a = 1.0
+--     b = 3.5
+--     eps = 5.0e-6
+--     u = x !! 0
+--     v = x !! 1
+--     w = x !! 2
+--
+-- lSaxis :: [[Double]] -> P.Axis B D.V2 Double
+-- lSaxis xs = P.r2Axis &~ do
+--   let ts = xs!!0
+--       us = xs!!1
+--       vs = xs!!2
+--       ws = xs!!3
+--   P.linePlot' $ zip ts us
+--   P.linePlot' $ zip ts vs
+--   P.linePlot' $ zip ts ws
+--
+-- main = do
+--   let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
+--   renderRasterific "diagrams/brusselator.png"
+--                    (D.dims2D 500.0 500.0)
+--                    (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
+-- @
+--
+-- With Sundials ARKode, it is possible to retrieve the Butcher tableau for the solver.
+--
+-- @
+-- import           Numeric.Sundials.ARKode.ODE
+-- import           Numeric.LinearAlgebra
+--
+-- import           Data.List (intercalate)
+--
+-- import           Text.PrettyPrint.HughesPJClass
+--
+--
+-- butcherTableauTex :: ButcherTable -> String
+-- butcherTableauTex (ButcherTable m c b b2) =
+--   render $
+--   vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}")
+--        , us
+--        , text "\\hline"
+--        , text bs <+> text "\\\\"
+--        , text b2s <+> text "\\\\"
+--        , text "\\end{array}"
+--        ]
+--   where
+--     n = rows m
+--     rs = toLists m
+--     ss = map (\r -> intercalate " & " $ map show r) rs
+--     ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss
+--     us = vcat $ map (\r -> text r <+> text "\\\\") ts
+--     bs  = " & " ++ (intercalate " & " $ map show $ toList b)
+--     b2s = " & " ++ (intercalate " & " $ map show $ toList b2)
+--
+-- main :: IO ()
+-- main = do
+--
+--   let res = butcherTable (SDIRK_2_1_2 undefined)
+--   putStrLn $ show res
+--   putStrLn $ butcherTableauTex res
+--
+--   let resA = butcherTable (KVAERNO_4_2_3 undefined)
+--   putStrLn $ show resA
+--   putStrLn $ butcherTableauTex resA
+--
+--   let resB = butcherTable (SDIRK_5_3_4 undefined)
+--   putStrLn $ show resB
+--   putStrLn $ butcherTableauTex resB
+-- @
+--
+--  Using the code above from the examples gives
+--
+-- KVAERNO_4_2_3
+--
+-- \[
+-- \begin{array}{c|cccc}
+-- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
+-- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\
+-- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
+-- 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
+-- \hline
+--  & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
+--  & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
+-- \end{array}
+-- \]
+--
+-- SDIRK_2_1_2
+--
+-- \[
+-- \begin{array}{c|cc}
+-- 1.0 & 1.0 & 0.0 \\
+-- 0.0 & -1.0 & 1.0 \\
+-- \hline
+--  & 0.5 & 0.5 \\
+--  & 1.0 & 0.0 \\
+-- \end{array}
+-- \]
+--
+-- SDIRK_5_3_4
+--
+-- \[
+-- \begin{array}{c|ccccc}
+-- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\
+-- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\
+-- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\
+-- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\
+-- 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
+-- \hline
+--  & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
+--  & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\
+-- \end{array}
+-- \]
+-----------------------------------------------------------------------------
+module Numeric.Sundials.ARKode.ODE ( odeSolve
+                                   , odeSolveV
+                                   , odeSolveVWith
+                                   , odeSolveVWith'
+                                   , ButcherTable(..)
+                                   , butcherTable
+                                   , ODEMethod(..)
+                                   , StepControl(..)
+                                   ) where
+
+import qualified Language.C.Inline as C
+import qualified Language.C.Inline.Unsafe as CU
+
+import           Data.Monoid ((<>))
+import           Data.Maybe (isJust)
+
+import           Foreign.C.Types (CDouble, CInt, CLong)
+import           Foreign.Ptr (Ptr)
+import           Foreign.Storable (poke)
+
+import qualified Data.Vector.Storable as V
+
+import           Data.Coerce (coerce)
+import           System.IO.Unsafe (unsafePerformIO)
+import           GHC.Generics (C1, Constructor, (:+:)(..), D1, Rep, Generic, M1(..),
+                               from, conName)
+
+import           Numeric.LinearAlgebra.Devel (createVector)
+
+import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
+                                                cols, toLists, size, reshape,
+                                                subVector, subMatrix, (><))
+
+import           Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..))
+import qualified Numeric.Sundials.Arkode as T
+import           Numeric.Sundials.Arkode (getDataFromContents, putDataInContents, arkSMax,
+                                          sDIRK_2_1_2,
+                                          bILLINGTON_3_3_2,
+                                          tRBDF2_3_3_2,
+                                          kVAERNO_4_2_3,
+                                          aRK324L2SA_DIRK_4_2_3,
+                                          cASH_5_2_4,
+                                          cASH_5_3_4,
+                                          sDIRK_5_3_4,
+                                          kVAERNO_5_3_4,
+                                          aRK436L2SA_DIRK_6_3_4,
+                                          kVAERNO_7_4_5,
+                                          aRK548L2SA_DIRK_8_4_5,
+                                          hEUN_EULER_2_1_2,
+                                          bOGACKI_SHAMPINE_4_2_3,
+                                          aRK324L2SA_ERK_4_2_3,
+                                          zONNEVELD_5_3_4,
+                                          aRK436L2SA_ERK_6_3_4,
+                                          sAYFY_ABURUB_6_3_4,
+                                          cASH_KARP_6_4_5,
+                                          fEHLBERG_6_4_5,
+                                          dORMAND_PRINCE_7_4_5,
+                                          aRK548L2SA_ERK_8_4_5,
+                                          vERNER_8_5_6,
+                                          fEHLBERG_13_7_8)
+
+
+C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
+
+C.include "<stdlib.h>"
+C.include "<stdio.h>"
+C.include "<math.h>"
+C.include "<arkode/arkode.h>"                 -- prototypes for ARKODE fcts., consts.
+C.include "<nvector/nvector_serial.h>"        -- serial N_Vector types, fcts., macros
+C.include "<sunmatrix/sunmatrix_dense.h>"     -- access to dense SUNMatrix
+C.include "<sunlinsol/sunlinsol_dense.h>"     -- access to dense SUNLinearSolver
+C.include "<arkode/arkode_direct.h>"          -- access to ARKDls interface
+C.include "<sundials/sundials_types.h>"       -- definition of type realtype
+C.include "<sundials/sundials_math.h>"
+C.include "../../../helpers.h"
+C.include "Numeric/Sundials/Arkode_hsc.h"
+
+
+-- | Stepping functions
+data ODEMethod = SDIRK_2_1_2            Jacobian
+               | SDIRK_2_1_2'
+               | BILLINGTON_3_3_2       Jacobian
+               | BILLINGTON_3_3_2'
+               | TRBDF2_3_3_2           Jacobian
+               | TRBDF2_3_3_2'
+               | KVAERNO_4_2_3          Jacobian
+               | KVAERNO_4_2_3'
+               | ARK324L2SA_DIRK_4_2_3  Jacobian
+               | ARK324L2SA_DIRK_4_2_3'
+               | CASH_5_2_4             Jacobian
+               | CASH_5_2_4'
+               | CASH_5_3_4             Jacobian
+               | CASH_5_3_4'
+               | SDIRK_5_3_4            Jacobian
+               | SDIRK_5_3_4'
+               | KVAERNO_5_3_4          Jacobian
+               | KVAERNO_5_3_4'
+               | ARK436L2SA_DIRK_6_3_4  Jacobian
+               | ARK436L2SA_DIRK_6_3_4'
+               | KVAERNO_7_4_5          Jacobian
+               | KVAERNO_7_4_5'
+               | ARK548L2SA_DIRK_8_4_5  Jacobian
+               | ARK548L2SA_DIRK_8_4_5'
+               | HEUN_EULER_2_1_2         Jacobian
+               | HEUN_EULER_2_1_2'
+               | BOGACKI_SHAMPINE_4_2_3   Jacobian
+               | BOGACKI_SHAMPINE_4_2_3'
+               | ARK324L2SA_ERK_4_2_3     Jacobian
+               | ARK324L2SA_ERK_4_2_3'
+               | ZONNEVELD_5_3_4          Jacobian
+               | ZONNEVELD_5_3_4'
+               | ARK436L2SA_ERK_6_3_4     Jacobian
+               | ARK436L2SA_ERK_6_3_4'
+               | SAYFY_ABURUB_6_3_4       Jacobian
+               | SAYFY_ABURUB_6_3_4'
+               | CASH_KARP_6_4_5          Jacobian
+               | CASH_KARP_6_4_5'
+               | FEHLBERG_6_4_5         Jacobian
+               | FEHLBERG_6_4_5'
+               | DORMAND_PRINCE_7_4_5     Jacobian
+               | DORMAND_PRINCE_7_4_5'
+               | ARK548L2SA_ERK_8_4_5     Jacobian
+               | ARK548L2SA_ERK_8_4_5'
+               | VERNER_8_5_6            Jacobian
+               | VERNER_8_5_6'
+               | FEHLBERG_13_7_8         Jacobian
+               | FEHLBERG_13_7_8'
+  deriving Generic
+
+constrName :: (HasConstructor (Rep a), Generic a)=> a -> String
+constrName = genericConstrName . from
+
+class HasConstructor (f :: * -> *) where
+  genericConstrName :: f x -> String
+
+instance HasConstructor f => HasConstructor (D1 c f) where
+  genericConstrName (M1 x) = genericConstrName x
+
+instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where
+  genericConstrName (L1 l) = genericConstrName l
+  genericConstrName (R1 r) = genericConstrName r
+
+instance Constructor c => HasConstructor (C1 c f) where
+  genericConstrName x = conName x
+
+instance Show ODEMethod where
+  show x = constrName x
+
+-- FIXME: We can probably do better here with generics
+getMethod :: ODEMethod -> Int
+getMethod (SDIRK_2_1_2 _)            = sDIRK_2_1_2
+getMethod (SDIRK_2_1_2')             = sDIRK_2_1_2
+getMethod (BILLINGTON_3_3_2 _)       = bILLINGTON_3_3_2
+getMethod (BILLINGTON_3_3_2')        = bILLINGTON_3_3_2
+getMethod (TRBDF2_3_3_2 _)           = tRBDF2_3_3_2
+getMethod (TRBDF2_3_3_2')            = tRBDF2_3_3_2
+getMethod (KVAERNO_4_2_3  _)         = kVAERNO_4_2_3
+getMethod (KVAERNO_4_2_3')           = kVAERNO_4_2_3
+getMethod (ARK324L2SA_DIRK_4_2_3 _)  = aRK324L2SA_DIRK_4_2_3
+getMethod (ARK324L2SA_DIRK_4_2_3')   = aRK324L2SA_DIRK_4_2_3
+getMethod (CASH_5_2_4 _)             = cASH_5_2_4
+getMethod (CASH_5_2_4')              = cASH_5_2_4
+getMethod (CASH_5_3_4 _)             = cASH_5_3_4
+getMethod (CASH_5_3_4')              = cASH_5_3_4
+getMethod (SDIRK_5_3_4 _)            = sDIRK_5_3_4
+getMethod (SDIRK_5_3_4')             = sDIRK_5_3_4
+getMethod (KVAERNO_5_3_4 _)          = kVAERNO_5_3_4
+getMethod (KVAERNO_5_3_4')           = kVAERNO_5_3_4
+getMethod (ARK436L2SA_DIRK_6_3_4 _)  = aRK436L2SA_DIRK_6_3_4
+getMethod (ARK436L2SA_DIRK_6_3_4')   = aRK436L2SA_DIRK_6_3_4
+getMethod (KVAERNO_7_4_5 _)          = kVAERNO_7_4_5
+getMethod (KVAERNO_7_4_5')           = kVAERNO_7_4_5
+getMethod (ARK548L2SA_DIRK_8_4_5 _)  = aRK548L2SA_DIRK_8_4_5
+getMethod (ARK548L2SA_DIRK_8_4_5')   = aRK548L2SA_DIRK_8_4_5
+getMethod (HEUN_EULER_2_1_2 _)       = hEUN_EULER_2_1_2
+getMethod (HEUN_EULER_2_1_2')        = hEUN_EULER_2_1_2
+getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3
+getMethod (BOGACKI_SHAMPINE_4_2_3')  = bOGACKI_SHAMPINE_4_2_3
+getMethod (ARK324L2SA_ERK_4_2_3 _)   = aRK324L2SA_ERK_4_2_3
+getMethod (ARK324L2SA_ERK_4_2_3')    = aRK324L2SA_ERK_4_2_3
+getMethod (ZONNEVELD_5_3_4 _)        = zONNEVELD_5_3_4
+getMethod (ZONNEVELD_5_3_4')         = zONNEVELD_5_3_4
+getMethod (ARK436L2SA_ERK_6_3_4 _)   = aRK436L2SA_ERK_6_3_4
+getMethod (ARK436L2SA_ERK_6_3_4')    = aRK436L2SA_ERK_6_3_4
+getMethod (SAYFY_ABURUB_6_3_4 _)     = sAYFY_ABURUB_6_3_4
+getMethod (SAYFY_ABURUB_6_3_4')      = sAYFY_ABURUB_6_3_4
+getMethod (CASH_KARP_6_4_5 _)        = cASH_KARP_6_4_5
+getMethod (CASH_KARP_6_4_5')         = cASH_KARP_6_4_5
+getMethod (FEHLBERG_6_4_5 _)         = fEHLBERG_6_4_5
+getMethod (FEHLBERG_6_4_5' )         = fEHLBERG_6_4_5
+getMethod (DORMAND_PRINCE_7_4_5 _)   = dORMAND_PRINCE_7_4_5
+getMethod (DORMAND_PRINCE_7_4_5')    = dORMAND_PRINCE_7_4_5
+getMethod (ARK548L2SA_ERK_8_4_5 _)   = aRK548L2SA_ERK_8_4_5
+getMethod (ARK548L2SA_ERK_8_4_5')    = aRK548L2SA_ERK_8_4_5
+getMethod (VERNER_8_5_6 _)           = vERNER_8_5_6
+getMethod (VERNER_8_5_6')            = vERNER_8_5_6
+getMethod (FEHLBERG_13_7_8 _)        = fEHLBERG_13_7_8
+getMethod (FEHLBERG_13_7_8')         = fEHLBERG_13_7_8
+
+getJacobian :: ODEMethod -> Maybe Jacobian
+getJacobian (SDIRK_2_1_2 j)            = Just j
+getJacobian (BILLINGTON_3_3_2 j)       = Just j
+getJacobian (TRBDF2_3_3_2 j)           = Just j
+getJacobian (KVAERNO_4_2_3  j)         = Just j
+getJacobian (ARK324L2SA_DIRK_4_2_3 j)  = Just j
+getJacobian (CASH_5_2_4 j)             = Just j
+getJacobian (CASH_5_3_4 j)             = Just j
+getJacobian (SDIRK_5_3_4 j)            = Just j
+getJacobian (KVAERNO_5_3_4 j)          = Just j
+getJacobian (ARK436L2SA_DIRK_6_3_4 j)  = Just j
+getJacobian (KVAERNO_7_4_5 j)          = Just j
+getJacobian (ARK548L2SA_DIRK_8_4_5 j)  = Just j
+getJacobian (HEUN_EULER_2_1_2 j)       = Just j
+getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j
+getJacobian (ARK324L2SA_ERK_4_2_3 j)   = Just j
+getJacobian (ZONNEVELD_5_3_4 j)        = Just j
+getJacobian (ARK436L2SA_ERK_6_3_4 j)   = Just j
+getJacobian (SAYFY_ABURUB_6_3_4 j)     = Just j
+getJacobian (CASH_KARP_6_4_5 j)        = Just j
+getJacobian (FEHLBERG_6_4_5 j)         = Just j
+getJacobian (DORMAND_PRINCE_7_4_5 j)   = Just j
+getJacobian (ARK548L2SA_ERK_8_4_5 j)   = Just j
+getJacobian (VERNER_8_5_6 j)           = Just j
+getJacobian (FEHLBERG_13_7_8 j)        = Just j
+getJacobian _                          = Nothing
+
+-- | A version of 'odeSolveVWith' with reasonable default step control.
+odeSolveV
+    :: ODEMethod
+    -> Maybe Double      -- ^ initial step size - by default, ARKode
+                         -- estimates the initial step size to be the
+                         -- solution \(h\) of the equation
+                         -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
+                         -- \(\ddot{y}\) is an estimated value of the
+                         -- second derivative of the solution at \(t_0\)
+    -> Double            -- ^ absolute tolerance for the state vector
+    -> Double            -- ^ relative tolerance for the state vector
+    -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+    -> Vector Double     -- ^ initial conditions
+    -> Vector Double     -- ^ desired solution times
+    -> Matrix Double     -- ^ solution
+odeSolveV meth hi epsAbs epsRel f y0 ts =
+  odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts
+    where
+      g t x0 = coerce $ f t x0
+
+-- | A version of 'odeSolveV' with reasonable default parameters and
+-- system of equations defined using lists. FIXME: we should say
+-- something about the fact we could use the Jacobian but don't for
+-- compatibility with hmatrix-gsl.
+odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+         -> [Double]                         -- ^ initial conditions
+         -> Vector Double                    -- ^ desired solution times
+         -> Matrix Double                    -- ^ solution
+odeSolve f y0 ts =
+  -- FIXME: These tolerances are different from the ones in GSL
+  odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1)  Nothing g (V.fromList y0) (V.fromList $ toList ts)
+  where
+    g t x0 = V.fromList $ f t (V.toList x0)
+
+odeSolveVWith ::
+  ODEMethod
+  -> StepControl
+  -> Maybe Double -- ^ initial step size - by default, ARKode
+                  -- estimates the initial step size to be the
+                  -- solution \(h\) of the equation
+                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
+                  -- \(\ddot{y}\) is an estimated value of the second
+                  -- derivative of the solution at \(t_0\)
+  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+  -> V.Vector Double                     -- ^ Initial conditions
+  -> V.Vector Double                     -- ^ Desired solution times
+  -> Matrix Double                       -- ^ Error code or solution
+odeSolveVWith method control initStepSize f y0 tt =
+  case odeSolveVWith' opts method control initStepSize f y0 tt of
+    Left c        -> error $ show c -- FIXME
+    Right (v, _d) -> v
+  where
+    opts = ODEOpts { maxNumSteps = 10000
+                   , minStep     = 1.0e-12
+                   , relTol      = error "relTol"
+                   , absTols     = error "absTol"
+                   , initStep    = error "initStep"
+                   , maxFail     = 10
+                   }
+
+odeSolveVWith' ::
+  ODEOpts
+  -> ODEMethod
+  -> StepControl
+  -> Maybe Double -- ^ initial step size - by default, ARKode
+                  -- estimates the initial step size to be the
+                  -- solution \(h\) of the equation
+                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
+                  -- \(\ddot{y}\) is an estimated value of the second
+                  -- derivative of the solution at \(t_0\)
+  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+  -> V.Vector Double                     -- ^ Initial conditions
+  -> V.Vector Double                     -- ^ Desired solution times
+  -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
+odeSolveVWith' opts method control initStepSize f y0 tt =
+  case solveOdeC (fromIntegral $ maxFail opts)
+                 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
+                 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
+                 (coerce f) (coerce y0) (coerce tt) of
+    Left c -> Left $ fromIntegral c
+    Right (v, d) ->  Right (reshape l (coerce v), d)
+  where
+    l = size y0
+    scise (X aTol rTol)                          = coerce (V.replicate l aTol, rTol)
+    scise (X' aTol rTol)                         = coerce (V.replicate l aTol, rTol)
+    scise (XX' aTol rTol yScale _yDotScale)      = coerce (V.replicate l aTol, yScale * rTol)
+    -- FIXME; Should we check that the length of ss is correct?
+    scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
+    jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
+           getJacobian method
+    matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
+      where
+        nr = fromIntegral $ rows m
+        nc = fromIntegral $ cols m
+        -- FIXME: efficiency
+        vs = V.fromList $ map coerce $ concat $ toLists m
+
+solveOdeC ::
+  CInt ->
+  CLong ->
+  CDouble ->
+  CInt ->
+  Maybe CDouble ->
+  (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
+  (V.Vector CDouble, CDouble) ->
+  (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+  -> V.Vector CDouble -- ^ Initial conditions
+  -> V.Vector CDouble -- ^ Desired solution times
+  -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
+solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
+          jacH (aTols, rTol) fun f0 ts = unsafePerformIO $ do
+
+  let isInitStepSize :: CInt
+      isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
+      ss :: CDouble
+      ss = case initStepSize of
+             -- It would be better to put an error message here but
+             -- inline-c seems to evaluate this even if it is never
+             -- used :(
+             Nothing -> 0.0
+             Just x  -> x
+
+  let dim = V.length f0
+      nEq :: CLong
+      nEq = fromIntegral dim
+      nTs :: CInt
+      nTs = fromIntegral $ V.length ts
+  -- FIXME: I believe this gets taken from the ghc heap and so should
+  -- be subject to garbage collection.
+  quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
+  qMatMut <- V.thaw quasiMatrixRes
+  diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
+  diagMut <- V.thaw diagnostics
+  -- We need the types that sundials expects. These are tied together
+  -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
+  let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
+      funIO x y f _ptr = do
+        -- Convert the pointer we get from C (y) to a vector, and then
+        -- apply the user-supplied function.
+        fImm <- fun x <$> getDataFromContents dim y
+        -- Fill in the provided pointer with the resulting vector.
+        putDataInContents fImm dim f
+        -- FIXME: I don't understand what this comment means
+        -- Unsafe since the function will be called many times.
+        [CU.exp| int{ 0 } |]
+  let isJac :: CInt
+      isJac = fromIntegral $ fromEnum $ isJust jacH
+      jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
+               Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
+               IO CInt
+      jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
+        case jacH of
+          Nothing   -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
+          Just jacI -> do j <- jacI t <$> getDataFromContents dim y
+                          poke jacS j
+                          -- FIXME: I don't understand what this comment means
+                          -- Unsafe since the function will be called many times.
+                          [CU.exp| int{ 0 } |]
+
+  res <- [C.block| int {
+                         /* general problem variables */
+
+                         int flag;                  /* reusable error-checking flag                 */
+                         int i, j;                  /* reusable loop indices                        */
+                         N_Vector y = NULL;         /* empty vector for storing solution            */
+                         N_Vector tv = NULL;        /* empty vector for storing absolute tolerances */
+                         SUNMatrix A = NULL;        /* empty matrix for linear solver               */
+                         SUNLinearSolver LS = NULL; /* empty linear solver object                   */
+                         void *arkode_mem = NULL;   /* empty ARKode memory structure                */
+                         realtype t;
+                         long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
+
+                         /* general problem parameters */
+
+                         realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time              */
+                         sunindextype NEQ = $(sunindextype nEq);             /* number of dependent vars. */
+
+                         /* Initialize data structures */
+
+                         y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
+                         if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
+                         /* Specify initial condition */
+                         for (i = 0; i < NEQ; i++) {
+                           NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
+                         };
+
+                         tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
+                         if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
+                         /* Specify tolerances */
+                         for (i = 0; i < NEQ; i++) {
+                           NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
+                         };
+
+                         arkode_mem = ARKodeCreate(); /* Create the solver memory */
+                         if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
+
+                         /* Call ARKodeInit to initialize the integrator memory and specify the */
+                         /* right-hand side function in y'=f(t,y), the inital time T0, and      */
+                         /* the initial dependent variable vector y.  Note: we treat the        */
+                         /* problem as fully implicit and set f_E to NULL and f_I to f.         */
+
+                         /* Here we use the C types defined in helpers.h which tie up with */
+                         /* the Haskell types defined in CLangToHaskellTypes                             */
+                         if ($(int method) < MIN_DIRK_NUM) {
+                           flag = ARKodeInit(arkode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y);
+                           if (check_flag(&flag, "ARKodeInit", 1)) return 1;
+                         } else {
+                           flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
+                           if (check_flag(&flag, "ARKodeInit", 1)) return 1;
+                         }
+
+                         flag = ARKodeSetMinStep(arkode_mem, $(double minStep_));
+                         if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1;
+                         flag = ARKodeSetMaxNumSteps(arkode_mem, $(long int maxNumSteps_));
+                         if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1;
+                         flag = ARKodeSetMaxErrTestFails(arkode_mem, $(int maxErrTestFails));
+                         if (check_flag(&flag, "ARKodeSetMaxErrTestFails", 1)) return 1;
+
+                         /* Set routines */
+                         flag = ARKodeSVtolerances(arkode_mem, $(double rTol), tv);
+                         if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1;
+
+                         /* Initialize dense matrix data structure and solver */
+                         A = SUNDenseMatrix(NEQ, NEQ);
+                         if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
+                         LS = SUNDenseLinearSolver(y, A);
+                         if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
+
+                         /* Attach matrix and linear solver */
+                         flag = ARKDlsSetLinearSolver(arkode_mem, LS, A);
+                         if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
+
+                         /* Set the initial step size if there is one */
+                         if ($(int isInitStepSize)) {
+                           /* FIXME: We could check if the initial step size is 0 */
+                           /* or even NaN and then throw an error                 */
+                           flag = ARKodeSetInitStep(arkode_mem, $(double ss));
+                           if (check_flag(&flag, "ARKodeSetInitStep", 1)) return 1;
+                         }
+
+                         /* Set the Jacobian if there is one */
+                         if ($(int isJac)) {
+                           flag = ARKDlsSetJacFn(arkode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
+                           if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
+                         }
+
+                         /* Store initial conditions */
+                         for (j = 0; j < NEQ; j++) {
+                           ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
+                         }
+
+                         /* Explicitly set the method */
+                         if ($(int method) >= MIN_DIRK_NUM) {
+                           flag = ARKodeSetIRKTableNum(arkode_mem, $(int method));
+                           if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
+                         } else {
+                           flag = ARKodeSetERKTableNum(arkode_mem, $(int method));
+                           if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
+                         }
+
+                         /* Main time-stepping loop: calls ARKode to perform the integration */
+                         /* Stops when the final time has been reached                       */
+                         for (i = 1; i < $(int nTs); i++) {
+
+                           flag = ARKode(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */
+                           if (check_flag(&flag, "ARKode", 1)) break;
+
+                           /* Store the results for Haskell */
+                           for (j = 0; j < NEQ; j++) {
+                             ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
+                           }
+
+                           /* unsuccessful solve: break */
+                           if (flag < 0) {
+                             fprintf(stderr,"Solver failure, stopping integration\n");
+                             break;
+                           }
+                         }
+
+                         /* Get some final statistics on how the solve progressed */
+
+                         flag = ARKodeGetNumSteps(arkode_mem, &nst);
+                         check_flag(&flag, "ARKodeGetNumSteps", 1);
+                         ($vec-ptr:(long int *diagMut))[0] = nst;
+
+                         flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
+                         check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
+                         ($vec-ptr:(long int *diagMut))[1] = nst_a;
+
+                         flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
+                         check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
+                         ($vec-ptr:(long int *diagMut))[2] = nfe;
+                         ($vec-ptr:(long int *diagMut))[3] = nfi;
+
+                         flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
+                         check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
+                         ($vec-ptr:(long int *diagMut))[4] = nsetups;
+
+                         flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
+                         check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
+                         ($vec-ptr:(long int *diagMut))[5] = netf;
+
+                         flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
+                         check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
+                         ($vec-ptr:(long int *diagMut))[6] = nni;
+
+                         flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
+                         check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
+                         ($vec-ptr:(long int *diagMut))[7] = ncfn;
+
+                         flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
+                         check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
+                         ($vec-ptr:(long int *diagMut))[8] = ncfn;
+
+                         flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
+                         check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
+                         ($vec-ptr:(long int *diagMut))[9] = ncfn;
+
+                         /* Clean up and return */
+                         N_VDestroy(y);            /* Free y vector          */
+                         N_VDestroy(tv);           /* Free tv vector         */
+                         ARKodeFree(&arkode_mem);  /* Free integrator memory */
+                         SUNLinSolFree(LS);        /* Free linear solver     */
+                         SUNMatDestroy(A);         /* Free A matrix          */
+
+                         return flag;
+                       } |]
+  if res == 0
+    then do
+      preD <- V.freeze diagMut
+      let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
+                                  (fromIntegral $ preD V.!1)
+                                  (fromIntegral $ preD V.!2)
+                                  (fromIntegral $ preD V.!3)
+                                  (fromIntegral $ preD V.!4)
+                                  (fromIntegral $ preD V.!5)
+                                  (fromIntegral $ preD V.!6)
+                                  (fromIntegral $ preD V.!7)
+                                  (fromIntegral $ preD V.!8)
+                                  (fromIntegral $ preD V.!9)
+      m <- V.freeze qMatMut
+      return $ Right (m, d)
+    else do
+      return $ Left res
+
+data ButcherTable = ButcherTable { am  :: Matrix Double
+                                 , cv  :: Vector Double
+                                 , bv  :: Vector Double
+                                 , b2v :: Vector Double
+                                 }
+  deriving Show
+
+data ButcherTable' a = ButcherTable' { am'  :: V.Vector a
+                                     , cv'  :: V.Vector a
+                                     , bv'  :: V.Vector a
+                                     , b2v' :: V.Vector a
+                                     }
+  deriving Show
+
+butcherTable :: ODEMethod -> ButcherTable
+butcherTable method =
+  case getBT method of
+    Left c -> error $ show c -- FIXME
+    Right (ButcherTable' v w x y, sqp) ->
+      ButcherTable { am = subMatrix (0, 0) (s, s) $ (arkSMax >< arkSMax) (V.toList v)
+                   , cv = subVector 0 s w
+                   , bv = subVector 0 s x
+                   , b2v = subVector 0 s y
+                   }
+      where
+        s = fromIntegral $ sqp V.! 0
+
+getBT :: ODEMethod -> Either Int (ButcherTable' Double, V.Vector Int)
+getBT method = case getButcherTable method of
+                 Left c ->
+                   Left $ fromIntegral c
+                 Right (ButcherTable' a b c d, sqp) ->
+                   Right $ ( ButcherTable' (coerce a) (coerce b) (coerce c) (coerce d)
+                           , V.map fromIntegral sqp )
+
+getButcherTable :: ODEMethod
+                -> Either CInt (ButcherTable' CDouble, V.Vector CInt)
+getButcherTable method = unsafePerformIO $ do
+  -- ARKode seems to want an ODE in order to set and then get the
+  -- Butcher tableau so here's one to keep it happy
+  let funI :: CDouble -> V.Vector CDouble -> V.Vector CDouble
+      funI _t ys = V.fromList [ ys V.! 0 ]
+  let funE :: CDouble -> V.Vector CDouble -> V.Vector CDouble
+      funE _t ys = V.fromList [ ys V.! 0 ]
+      f0        = V.fromList [ 1.0 ]
+      ts        = V.fromList [ 0.0 ]
+      dim = V.length f0
+      nEq :: CLong
+      nEq = fromIntegral dim
+      mN :: CInt
+      mN = fromIntegral $ getMethod method
+
+  btSQP :: V.Vector CInt <- createVector 3
+  btSQPMut <- V.thaw btSQP
+  btAs :: V.Vector CDouble <- createVector (arkSMax * arkSMax)
+  btAsMut <- V.thaw btAs
+  btCs  :: V.Vector CDouble <- createVector arkSMax
+  btBs  :: V.Vector CDouble <- createVector arkSMax
+  btB2s :: V.Vector CDouble <- createVector arkSMax
+  btCsMut  <- V.thaw btCs
+  btBsMut  <- V.thaw btBs
+  btB2sMut <- V.thaw btB2s
+  let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
+      funIOI x y f _ptr = do
+        fImm <- funI x <$> getDataFromContents dim y
+        putDataInContents fImm dim f
+        -- FIXME: I don't understand what this comment means
+        -- Unsafe since the function will be called many times.
+        [CU.exp| int{ 0 } |]
+  let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
+      funIOE x y f _ptr = do
+        fImm <- funE x <$> getDataFromContents dim y
+        putDataInContents fImm dim f
+        -- FIXME: I don't understand what this comment means
+        -- Unsafe since the function will be called many times.
+        [CU.exp| int{ 0 } |]
+  res <- [C.block| int {
+                         /* general problem variables */
+
+                         int flag;                /* reusable error-checking flag      */
+                         N_Vector y = NULL;       /* empty vector for storing solution */
+                         void *arkode_mem = NULL; /* empty ARKode memory structure     */
+                         int i, j;                /* reusable loop indices             */
+
+                         /* general problem parameters */
+
+                         realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time             */
+                         sunindextype NEQ = $(sunindextype nEq);           /* number of dependent vars */
+
+                         /* Initialize data structures */
+
+                         y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
+                         if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
+                         /* Specify initial condition */
+                         for (i = 0; i < NEQ; i++) {
+                           NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
+                         };
+                         arkode_mem = ARKodeCreate(); /* Create the solver memory */
+                         if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
+
+                         flag = ARKodeInit(arkode_mem, $fun:(int (* funIOE) (double t, SunVector y[], SunVector dydt[], void * params)), $fun:(int (* funIOI) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
+                         if (check_flag(&flag, "ARKodeInit", 1)) return 1;
+
+                         if ($(int mN) >= MIN_DIRK_NUM) {
+                         flag = ARKodeSetIRKTableNum(arkode_mem, $(int mN));
+                         if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
+                         } else {
+                         flag = ARKodeSetERKTableNum(arkode_mem, $(int mN));
+                         if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
+                         }
+
+                         int s, q, p;
+                         realtype *ai = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
+                         realtype *ae = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
+                         realtype *ci = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *ce = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *bi = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *be = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *b2i = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         realtype *b2e = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
+                         flag = ARKodeGetCurrentButcherTables(arkode_mem, &s, &q, &p, ai, ae, ci, ce, bi, be, b2i, b2e);
+                         if (check_flag(&flag, "ARKode", 1)) return 1;
+                         $vec-ptr:(int *btSQPMut)[0] = s;
+                         $vec-ptr:(int *btSQPMut)[1] = q;
+                         $vec-ptr:(int *btSQPMut)[2] = p;
+                         for (i = 0; i < s; i++) {
+                           for (j = 0; j < s; j++) {
+                             /* FIXME: double should be realtype */
+                             ($vec-ptr:(double *btAsMut))[i * ARK_S_MAX + j] = ai[i * ARK_S_MAX + j];
+                           }
+                         }
+
+                         for (i = 0; i < s; i++) {
+                           ($vec-ptr:(double *btCsMut))[i]  = ci[i];
+                           ($vec-ptr:(double *btBsMut))[i]  = bi[i];
+                           ($vec-ptr:(double *btB2sMut))[i] = b2i[i];
+                         }
+
+                         /* Clean up and return */
+                         N_VDestroy(y);            /* Free y vector */
+                         ARKodeFree(&arkode_mem);  /* Free integrator memory */
+
+                         return flag;
+                       } |]
+  if res == 0
+    then do
+      x <- V.freeze btAsMut
+      y <- V.freeze btSQPMut
+      z <- V.freeze btCsMut
+      u <- V.freeze btBsMut
+      v <- V.freeze btB2sMut
+      return $ Right (ButcherTable' { am' = x, cv' = z, bv' = u, b2v' = v }, y)
+    else do
+      return $ Left res
+
+-- | Adaptive step-size control
+-- functions.
+--
+-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
+-- allows the user to control the step size adjustment using
+-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
+-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
+-- is the required relative error, \(s_i\) is a vector of scaling
+-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
+-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
+--
+-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
+-- allows the user to control the step size adjustment using
+-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
+-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
+-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
+-- effect.
+data StepControl = X     Double Double -- ^ absolute and relative tolerance for \(y\); in GSL terms, \(a_{y} = 1\) and \(a_{dy/dt} = 0\); in ARKode terms, the \(\eta^{abs}_i\) are identical
+                 | X'    Double Double -- ^ absolute and relative tolerance for \(\dot{y}\); in GSL terms, \(a_{y} = 0\) and \(a_{dy/dt} = 1\); in ARKode terms, the latter is treated as the relative tolerance for \(y\) so this is the same as specifying 'X' which may be entirely incorrect for the given problem
+                 | XX'   Double Double Double Double -- ^ include both via relative tolerance
+                                                     -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
+                 | ScXX' Double Double Double Double (Vector Double) -- ^ scale absolute tolerance of \(y_i\); in ARKode terms, \(a_{{dy}/{dt}}\) is ignored, \(\eta^{abs}_i = s_i \epsilon^{abs}\) and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
diff --git a/packages/sundials/src/Numeric/Sundials/Arkode.hsc b/packages/sundials/src/Numeric/Sundials/Arkode.hsc
new file mode 100644
index 0000000..0850258
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/Arkode.hsc
@@ -0,0 +1,204 @@
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE EmptyDataDecls #-}
+
+module Numeric.Sundials.Arkode where
+
+import           Foreign
+import           Foreign.C.Types
+
+import           Language.C.Types as CT
+
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Storable.Mutable as VM
+
+import qualified Language.Haskell.TH as TH
+import qualified Data.Map as Map
+import           Language.C.Inline.Context
+
+import qualified Data.Vector.Storable as V
+
+
+#include <stdio.h>
+#include <sundials/sundials_nvector.h>
+#include <sundials/sundials_matrix.h>
+#include <nvector/nvector_serial.h>
+#include <sunmatrix/sunmatrix_dense.h>
+#include <arkode/arkode.h>
+#include <cvode/cvode.h>
+
+
+data SunVector
+data SunMatrix = SunMatrix { rows :: CInt
+                           , cols :: CInt
+                           , vals :: V.Vector CDouble
+                           }
+
+-- | This is true only if configured/ built as 64 bits
+type SunIndexType = CLong
+
+sunTypesTable :: Map.Map TypeSpecifier TH.TypeQ
+sunTypesTable = Map.fromList
+  [
+    (TypeName "sunindextype", [t| SunIndexType |] )
+  , (TypeName "SunVector",    [t| SunVector |] )
+  , (TypeName "SunMatrix",    [t| SunMatrix |] )
+  ]
+
+sunCtx :: Context
+sunCtx = mempty {ctxTypesTable = sunTypesTable}
+
+getMatrixDataFromContents :: Ptr SunMatrix -> IO SunMatrix
+getMatrixDataFromContents ptr = do
+  qtr <- getContentMatrixPtr ptr
+  rs  <- getNRows qtr
+  cs  <- getNCols qtr
+  rtr <- getMatrixData qtr
+  vs  <- vectorFromC (fromIntegral $ rs * cs) rtr
+  return $ SunMatrix { rows = rs, cols = cs, vals = vs }
+
+putMatrixDataFromContents :: SunMatrix -> Ptr SunMatrix -> IO ()
+putMatrixDataFromContents mat ptr = do
+  let rs = rows mat
+      cs = cols mat
+      vs = vals mat
+  qtr <- getContentMatrixPtr ptr
+  putNRows rs qtr
+  putNCols cs qtr
+  rtr <- getMatrixData qtr
+  vectorToC vs (fromIntegral $ rs * cs) rtr
+
+instance Storable SunMatrix where
+  poke        = flip putMatrixDataFromContents
+  peek        = getMatrixDataFromContents
+  sizeOf _    = error "sizeOf not supported for SunMatrix"
+  alignment _ = error "alignment not supported for SunMatrix"
+
+vectorFromC :: Storable a => Int -> Ptr a -> IO (VS.Vector a)
+vectorFromC len ptr = do
+  ptr' <- newForeignPtr_ ptr
+  VS.freeze $ VM.unsafeFromForeignPtr0 ptr' len
+
+vectorToC :: Storable a => VS.Vector a -> Int -> Ptr a -> IO ()
+vectorToC vec len ptr = do
+  ptr' <- newForeignPtr_ ptr
+  VS.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
+
+getDataFromContents :: Int -> Ptr SunVector -> IO (VS.Vector CDouble)
+getDataFromContents len ptr = do
+  qtr <- getContentPtr ptr
+  rtr <- getData qtr
+  vectorFromC len rtr
+
+putDataInContents :: Storable a => VS.Vector a -> Int -> Ptr b -> IO ()
+putDataInContents vec len ptr = do
+  qtr <- getContentPtr ptr
+  rtr <- getData qtr
+  vectorToC vec len rtr
+
+#def typedef struct _generic_N_Vector SunVector;
+#def typedef struct _N_VectorContent_Serial SunContent;
+
+#def typedef struct _generic_SUNMatrix SunMatrix;
+#def typedef struct _SUNMatrixContent_Dense SunMatrixContent;
+
+getContentMatrixPtr :: Storable a => Ptr b -> IO a
+getContentMatrixPtr ptr = (#peek SunMatrix, content) ptr
+
+getNRows :: Ptr b -> IO CInt
+getNRows ptr = (#peek SunMatrixContent, M) ptr
+putNRows :: CInt -> Ptr b -> IO ()
+putNRows nr ptr = (#poke SunMatrixContent, M) ptr nr
+
+getNCols :: Ptr b -> IO CInt
+getNCols ptr = (#peek SunMatrixContent, N) ptr
+putNCols :: CInt -> Ptr b -> IO ()
+putNCols nc ptr = (#poke SunMatrixContent, N) ptr nc
+
+getMatrixData :: Storable a => Ptr b -> IO a
+getMatrixData ptr = (#peek SunMatrixContent, data) ptr
+
+getContentPtr :: Storable a => Ptr b -> IO a
+getContentPtr ptr = (#peek SunVector, content) ptr
+
+getData :: Storable a => Ptr b -> IO a
+getData ptr = (#peek SunContent, data) ptr
+
+cV_ADAMS :: Int
+cV_ADAMS = #const CV_ADAMS
+cV_BDF :: Int
+cV_BDF = #const CV_BDF
+
+arkSMax :: Int
+arkSMax = #const ARK_S_MAX
+
+mIN_DIRK_NUM, mAX_DIRK_NUM :: Int
+mIN_DIRK_NUM = #const MIN_DIRK_NUM
+mAX_DIRK_NUM = #const MAX_DIRK_NUM
+
+-- FIXME: We could just use inline-c instead
+
+-- Butcher table accessors -- implicit
+sDIRK_2_1_2 :: Int
+sDIRK_2_1_2 = #const SDIRK_2_1_2
+bILLINGTON_3_3_2 :: Int
+bILLINGTON_3_3_2 = #const BILLINGTON_3_3_2
+tRBDF2_3_3_2 :: Int
+tRBDF2_3_3_2 = #const TRBDF2_3_3_2
+kVAERNO_4_2_3 :: Int
+kVAERNO_4_2_3 = #const KVAERNO_4_2_3
+aRK324L2SA_DIRK_4_2_3 :: Int
+aRK324L2SA_DIRK_4_2_3 = #const ARK324L2SA_DIRK_4_2_3
+cASH_5_2_4 :: Int
+cASH_5_2_4 = #const CASH_5_2_4
+cASH_5_3_4 :: Int
+cASH_5_3_4 = #const CASH_5_3_4
+sDIRK_5_3_4 :: Int
+sDIRK_5_3_4 = #const SDIRK_5_3_4
+kVAERNO_5_3_4 :: Int
+kVAERNO_5_3_4 = #const KVAERNO_5_3_4
+aRK436L2SA_DIRK_6_3_4 :: Int
+aRK436L2SA_DIRK_6_3_4 = #const ARK436L2SA_DIRK_6_3_4
+kVAERNO_7_4_5 :: Int
+kVAERNO_7_4_5 = #const KVAERNO_7_4_5
+aRK548L2SA_DIRK_8_4_5 :: Int
+aRK548L2SA_DIRK_8_4_5 = #const ARK548L2SA_DIRK_8_4_5
+
+-- #define DEFAULT_DIRK_2          SDIRK_2_1_2
+-- #define DEFAULT_DIRK_3          ARK324L2SA_DIRK_4_2_3
+-- #define DEFAULT_DIRK_4          SDIRK_5_3_4
+-- #define DEFAULT_DIRK_5          ARK548L2SA_DIRK_8_4_5
+
+-- Butcher table accessors -- explicit
+hEUN_EULER_2_1_2 :: Int
+hEUN_EULER_2_1_2 = #const HEUN_EULER_2_1_2
+bOGACKI_SHAMPINE_4_2_3 :: Int
+bOGACKI_SHAMPINE_4_2_3 = #const BOGACKI_SHAMPINE_4_2_3
+aRK324L2SA_ERK_4_2_3 :: Int
+aRK324L2SA_ERK_4_2_3 = #const ARK324L2SA_ERK_4_2_3
+zONNEVELD_5_3_4 :: Int
+zONNEVELD_5_3_4 = #const ZONNEVELD_5_3_4
+aRK436L2SA_ERK_6_3_4 :: Int
+aRK436L2SA_ERK_6_3_4 = #const ARK436L2SA_ERK_6_3_4
+sAYFY_ABURUB_6_3_4 :: Int
+sAYFY_ABURUB_6_3_4 = #const SAYFY_ABURUB_6_3_4
+cASH_KARP_6_4_5 :: Int
+cASH_KARP_6_4_5 = #const CASH_KARP_6_4_5
+fEHLBERG_6_4_5 :: Int
+fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5
+dORMAND_PRINCE_7_4_5 :: Int
+dORMAND_PRINCE_7_4_5 = #const DORMAND_PRINCE_7_4_5
+aRK548L2SA_ERK_8_4_5 :: Int
+aRK548L2SA_ERK_8_4_5 = #const ARK548L2SA_ERK_8_4_5
+vERNER_8_5_6 :: Int
+vERNER_8_5_6 = #const VERNER_8_5_6
+fEHLBERG_13_7_8 :: Int
+fEHLBERG_13_7_8 = #const FEHLBERG_13_7_8
+
+-- #define DEFAULT_ERK_2           HEUN_EULER_2_1_2
+-- #define DEFAULT_ERK_3           BOGACKI_SHAMPINE_4_2_3
+-- #define DEFAULT_ERK_4           ZONNEVELD_5_3_4
+-- #define DEFAULT_ERK_5           CASH_KARP_6_4_5
+-- #define DEFAULT_ERK_6           VERNER_8_5_6
+-- #define DEFAULT_ERK_8           FEHLBERG_13_7_8
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
new file mode 100644
index 0000000..a6f185e
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
@@ -0,0 +1,476 @@
+{-# OPTIONS_GHC -Wall #-}
+
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE MultiWayIf #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Sundials.CVode.ODE
+-- Copyright   :  Dominic Steinitz 2018,
+--                Novadiscovery 2018
+-- License     :  BSD
+-- Maintainer  :  Dominic Steinitz
+-- Stability   :  provisional
+--
+-- Solution of ordinary differential equation (ODE) initial value problems.
+--
+-- <https://computation.llnl.gov/projects/sundials/sundials-software>
+--
+-- A simple example:
+--
+-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
+--
+-- @
+-- import           Numeric.Sundials.CVode.ODE
+-- import           Numeric.LinearAlgebra
+--
+-- import           Plots as P
+-- import qualified Diagrams.Prelude as D
+-- import           Diagrams.Backend.Rasterific
+--
+-- brusselator :: Double -> [Double] -> [Double]
+-- brusselator _t x = [ a - (w + 1) * u + v * u * u
+--                    , w * u - v * u * u
+--                    , (b - w) / eps - w * u
+--                    ]
+--   where
+--     a = 1.0
+--     b = 3.5
+--     eps = 5.0e-6
+--     u = x !! 0
+--     v = x !! 1
+--     w = x !! 2
+--
+-- lSaxis :: [[Double]] -> P.Axis B D.V2 Double
+-- lSaxis xs = P.r2Axis &~ do
+--   let ts = xs!!0
+--       us = xs!!1
+--       vs = xs!!2
+--       ws = xs!!3
+--   P.linePlot' $ zip ts us
+--   P.linePlot' $ zip ts vs
+--   P.linePlot' $ zip ts ws
+--
+-- main = do
+--   let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
+--   renderRasterific "diagrams/brusselator.png"
+--                    (D.dims2D 500.0 500.0)
+--                    (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
+-- @
+--
+-----------------------------------------------------------------------------
+module Numeric.Sundials.CVode.ODE ( odeSolve
+                                   , odeSolveV
+                                   , odeSolveVWith
+                                   , odeSolveVWith'
+                                   , ODEMethod(..)
+                                   , StepControl(..)
+                                   ) where
+
+import qualified Language.C.Inline as C
+import qualified Language.C.Inline.Unsafe as CU
+
+import           Data.Monoid ((<>))
+import           Data.Maybe (isJust)
+
+import           Foreign.C.Types (CDouble, CInt, CLong)
+import           Foreign.Ptr (Ptr)
+import           Foreign.Storable (poke)
+
+import qualified Data.Vector.Storable as V
+
+import           Data.Coerce (coerce)
+import           System.IO.Unsafe (unsafePerformIO)
+
+import           Numeric.LinearAlgebra.Devel (createVector)
+
+import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
+                                                cols, toLists, size, reshape)
+
+import           Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF,
+                                          getDataFromContents, putDataInContents)
+import qualified Numeric.Sundials.Arkode as T
+import           Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..))
+
+
+C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
+
+C.include "<stdlib.h>"
+C.include "<stdio.h>"
+C.include "<math.h>"
+C.include "<cvode/cvode.h>"               -- prototypes for CVODE fcts., consts.
+C.include "<nvector/nvector_serial.h>"    -- serial N_Vector types, fcts., macros
+C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
+C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
+C.include "<cvode/cvode_direct.h>"        -- access to CVDls interface
+C.include "<sundials/sundials_types.h>"   -- definition of type realtype
+C.include "<sundials/sundials_math.h>"
+C.include "../../../helpers.h"
+C.include "Numeric/Sundials/Arkode_hsc.h"
+
+
+-- | Stepping functions
+data ODEMethod = ADAMS
+               | BDF
+
+getMethod :: ODEMethod -> Int
+getMethod (ADAMS) = cV_ADAMS
+getMethod (BDF)   = cV_BDF
+
+getJacobian :: ODEMethod -> Maybe Jacobian
+getJacobian _ = Nothing
+
+-- | A version of 'odeSolveVWith' with reasonable default step control.
+odeSolveV
+    :: ODEMethod
+    -> Maybe Double      -- ^ initial step size - by default, CVode
+                         -- estimates the initial step size to be the
+                         -- solution \(h\) of the equation
+                         -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
+                         -- \(\ddot{y}\) is an estimated value of the
+                         -- second derivative of the solution at \(t_0\)
+    -> Double            -- ^ absolute tolerance for the state vector
+    -> Double            -- ^ relative tolerance for the state vector
+    -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+    -> Vector Double     -- ^ initial conditions
+    -> Vector Double     -- ^ desired solution times
+    -> Matrix Double     -- ^ solution
+odeSolveV meth hi epsAbs epsRel f y0 ts =
+  odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts
+  where
+    g t x0 = coerce $ f t x0
+
+-- | A version of 'odeSolveV' with reasonable default parameters and
+-- system of equations defined using lists. FIXME: we should say
+-- something about the fact we could use the Jacobian but don't for
+-- compatibility with hmatrix-gsl.
+odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+         -> [Double]                         -- ^ initial conditions
+         -> Vector Double                    -- ^ desired solution times
+         -> Matrix Double                    -- ^ solution
+odeSolve f y0 ts =
+  -- FIXME: These tolerances are different from the ones in GSL
+  odeSolveVWith BDF (XX' 1.0e-6 1.0e-10 1 1)  Nothing g (V.fromList y0) (V.fromList $ toList ts)
+  where
+    g t x0 = V.fromList $ f t (V.toList x0)
+
+odeSolveVWith ::
+  ODEMethod
+  -> StepControl
+  -> Maybe Double -- ^ initial step size - by default, CVode
+                  -- estimates the initial step size to be the
+                  -- solution \(h\) of the equation
+                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
+                  -- \(\ddot{y}\) is an estimated value of the second
+                  -- derivative of the solution at \(t_0\)
+  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+  -> V.Vector Double                     -- ^ Initial conditions
+  -> V.Vector Double                     -- ^ Desired solution times
+  -> Matrix Double                       -- ^ Error code or solution
+odeSolveVWith method control initStepSize f y0 tt =
+  case odeSolveVWith' opts method control initStepSize f y0 tt of
+    Left c        -> error $ show c -- FIXME
+    Right (v, _d) -> v
+  where
+    opts = ODEOpts { maxNumSteps = 10000
+                   , minStep     = 1.0e-12
+                   , relTol      = error "relTol"
+                   , absTols     = error "absTol"
+                   , initStep    = error "initStep"
+                   , maxFail     = 10
+                   }
+
+odeSolveVWith' ::
+  ODEOpts
+  -> ODEMethod
+  -> StepControl
+  -> Maybe Double -- ^ initial step size - by default, CVode
+                  -- estimates the initial step size to be the
+                  -- solution \(h\) of the equation
+                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
+                  -- \(\ddot{y}\) is an estimated value of the second
+                  -- derivative of the solution at \(t_0\)
+  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+  -> V.Vector Double                     -- ^ Initial conditions
+  -> V.Vector Double                     -- ^ Desired solution times
+  -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
+odeSolveVWith' opts method control initStepSize f y0 tt =
+  case solveOdeC (fromIntegral $ maxFail opts)
+                 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
+                 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
+                 (coerce f) (coerce y0) (coerce tt) of
+    Left c -> Left $ fromIntegral c
+    Right (v, d) -> Right (reshape l (coerce v), d)
+  where
+    l = size y0
+    scise (X aTol rTol)                          = coerce (V.replicate l aTol, rTol)
+    scise (X' aTol rTol)                         = coerce (V.replicate l aTol, rTol)
+    scise (XX' aTol rTol yScale _yDotScale)      = coerce (V.replicate l aTol, yScale * rTol)
+    -- FIXME; Should we check that the length of ss is correct?
+    scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
+    jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
+           getJacobian method
+    matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
+      where
+        nr = fromIntegral $ rows m
+        nc = fromIntegral $ cols m
+        -- FIXME: efficiency
+        vs = V.fromList $ map coerce $ concat $ toLists m
+
+solveOdeC ::
+  CInt ->
+  CLong ->
+  CDouble ->
+  CInt ->
+  Maybe CDouble ->
+  (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
+  (V.Vector CDouble, CDouble) ->
+  (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
+  -> V.Vector CDouble -- ^ Initial conditions
+  -> V.Vector CDouble -- ^ Desired solution times
+  -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
+solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
+          jacH (aTols, rTol) fun f0 ts =
+  unsafePerformIO $ do
+
+  let isInitStepSize :: CInt
+      isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
+      ss :: CDouble
+      ss = case initStepSize of
+             -- It would be better to put an error message here but
+             -- inline-c seems to evaluate this even if it is never
+             -- used :(
+             Nothing -> 0.0
+             Just x  -> x
+
+  let dim = V.length f0
+      nEq :: CLong
+      nEq = fromIntegral dim
+      nTs :: CInt
+      nTs = fromIntegral $ V.length ts
+  quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
+  qMatMut <- V.thaw quasiMatrixRes
+  diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
+  diagMut <- V.thaw diagnostics
+  -- We need the types that sundials expects. These are tied together
+  -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
+  let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
+      funIO x y f _ptr = do
+        -- Convert the pointer we get from C (y) to a vector, and then
+        -- apply the user-supplied function.
+        fImm <- fun x <$> getDataFromContents dim y
+        -- Fill in the provided pointer with the resulting vector.
+        putDataInContents fImm dim f
+        -- FIXME: I don't understand what this comment means
+        -- Unsafe since the function will be called many times.
+        [CU.exp| int{ 0 } |]
+  let isJac :: CInt
+      isJac = fromIntegral $ fromEnum $ isJust jacH
+      jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
+               Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
+               IO CInt
+      jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
+        case jacH of
+          Nothing   -> error "Numeric.Sundials.CVode.ODE: Jacobian not defined"
+          Just jacI -> do j <- jacI t <$> getDataFromContents dim y
+                          poke jacS j
+                          -- FIXME: I don't understand what this comment means
+                          -- Unsafe since the function will be called many times.
+                          [CU.exp| int{ 0 } |]
+
+  res <- [C.block| int {
+                         /* general problem variables */
+
+                         int flag;                  /* reusable error-checking flag                 */
+                         int i, j;                  /* reusable loop indices                        */
+                         N_Vector y = NULL;         /* empty vector for storing solution            */
+                         N_Vector tv = NULL;        /* empty vector for storing absolute tolerances */
+
+                         SUNMatrix A = NULL;        /* empty matrix for linear solver               */
+                         SUNLinearSolver LS = NULL; /* empty linear solver object                   */
+                         void *cvode_mem = NULL;    /* empty CVODE memory structure                 */
+                         realtype t;
+                         long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
+
+                         /* general problem parameters */
+
+                         realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time              */
+                         sunindextype NEQ = $(sunindextype nEq);           /* number of dependent vars. */
+
+                         /* Initialize data structures */
+
+                         y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
+                         if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
+                         /* Specify initial condition */
+                         for (i = 0; i < NEQ; i++) {
+                           NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
+                         };
+
+                         cvode_mem = CVodeCreate($(int method), CV_NEWTON);
+                         if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
+
+                         /* Call CVodeInit to initialize the integrator memory and specify the
+                          * user's right hand side function in y'=f(t,y), the inital time T0, and
+                          * the initial dependent variable vector y. */
+                         flag = CVodeInit(cvode_mem,   $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
+                         if (check_flag(&flag, "CVodeInit", 1)) return(1);
+
+                         tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
+                         if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
+                         /* Specify tolerances */
+                         for (i = 0; i < NEQ; i++) {
+                           NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
+                         };
+
+                         flag = CVodeSetMinStep(cvode_mem, $(double minStep_));
+                         if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1;
+                         flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_));
+                         if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1;
+                         flag = CVodeSetMaxErrTestFails(cvode_mem, $(int maxErrTestFails));
+                         if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return 1;
+
+                         /* Call CVodeSVtolerances to specify the scalar relative tolerance
+                          * and vector absolute tolerances */
+                         flag = CVodeSVtolerances(cvode_mem, $(double rTol), tv);
+                         if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
+
+                         /* Initialize dense matrix data structure and solver */
+                         A = SUNDenseMatrix(NEQ, NEQ);
+                         if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
+                         LS = SUNDenseLinearSolver(y, A);
+                         if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
+
+                         /* Attach matrix and linear solver */
+                         flag = CVDlsSetLinearSolver(cvode_mem, LS, A);
+                         if (check_flag(&flag, "CVDlsSetLinearSolver", 1)) return 1;
+
+                         /* Set the initial step size if there is one */
+                         if ($(int isInitStepSize)) {
+                           /* FIXME: We could check if the initial step size is 0 */
+                           /* or even NaN and then throw an error                 */
+                           flag = CVodeSetInitStep(cvode_mem, $(double ss));
+                           if (check_flag(&flag, "CVodeSetInitStep", 1)) return 1;
+                         }
+
+                         /* Set the Jacobian if there is one */
+                         if ($(int isJac)) {
+                           flag = CVDlsSetJacFn(cvode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
+                           if (check_flag(&flag, "CVDlsSetJacFn", 1)) return 1;
+                         }
+
+                         /* Store initial conditions */
+                         for (j = 0; j < NEQ; j++) {
+                           ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
+                         }
+
+                         /* Main time-stepping loop: calls CVode to perform the integration */
+                         /* Stops when the final time has been reached                       */
+                         for (i = 1; i < $(int nTs); i++) {
+
+                           flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[i], y, &t, CV_NORMAL); /* call integrator */
+                           if (check_flag(&flag, "CVode", 1)) break;
+
+                           /* Store the results for Haskell */
+                           for (j = 0; j < NEQ; j++) {
+                             ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
+                           }
+
+                           /* unsuccessful solve: break */
+                           if (flag < 0) {
+                             fprintf(stderr,"Solver failure, stopping integration\n");
+                             break;
+                           }
+                         }
+
+                         /* Get some final statistics on how the solve progressed */
+
+                         flag = CVodeGetNumSteps(cvode_mem, &nst);
+                         check_flag(&flag, "CVodeGetNumSteps", 1);
+                         ($vec-ptr:(long int *diagMut))[0] = nst;
+
+                         /* FIXME */
+                         ($vec-ptr:(long int *diagMut))[1] = 0;
+
+                         flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
+                         check_flag(&flag, "CVodeGetNumRhsEvals", 1);
+                         ($vec-ptr:(long int *diagMut))[2] = nfe;
+                         /* FIXME */
+                         ($vec-ptr:(long int *diagMut))[3] = 0;
+
+                         flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
+                         check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
+                         ($vec-ptr:(long int *diagMut))[4] = nsetups;
+
+                         flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
+                         check_flag(&flag, "CVodeGetNumErrTestFails", 1);
+                         ($vec-ptr:(long int *diagMut))[5] = netf;
+
+                         flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
+                         check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
+                         ($vec-ptr:(long int *diagMut))[6] = nni;
+
+                         flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
+                         check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
+                         ($vec-ptr:(long int *diagMut))[7] = ncfn;
+
+                         flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
+                         check_flag(&flag, "CVDlsGetNumJacEvals", 1);
+                         ($vec-ptr:(long int *diagMut))[8] = ncfn;
+
+                         flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
+                         check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
+                         ($vec-ptr:(long int *diagMut))[9] = ncfn;
+
+                         /* Clean up and return */
+
+                         N_VDestroy(y);          /* Free y vector          */
+                         N_VDestroy(tv);         /* Free tv vector         */
+                         CVodeFree(&cvode_mem);  /* Free integrator memory */
+                         SUNLinSolFree(LS);      /* Free linear solver     */
+                         SUNMatDestroy(A);       /* Free A matrix          */
+
+                         return flag;
+                       } |]
+  if res == 0
+    then do
+      preD <- V.freeze diagMut
+      let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
+                                  (fromIntegral $ preD V.!1)
+                                  (fromIntegral $ preD V.!2)
+                                  (fromIntegral $ preD V.!3)
+                                  (fromIntegral $ preD V.!4)
+                                  (fromIntegral $ preD V.!5)
+                                  (fromIntegral $ preD V.!6)
+                                  (fromIntegral $ preD V.!7)
+                                  (fromIntegral $ preD V.!8)
+                                  (fromIntegral $ preD V.!9)
+      m <- V.freeze qMatMut
+      return $ Right (m, d)
+    else do
+      return $ Left res
+
+-- | Adaptive step-size control
+-- functions.
+--
+-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
+-- allows the user to control the step size adjustment using
+-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
+-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
+-- is the required relative error, \(s_i\) is a vector of scaling
+-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
+-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
+--
+-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
+-- allows the user to control the step size adjustment using
+-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
+-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
+-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
+-- effect.
+data StepControl = X     Double Double -- ^ absolute and relative tolerance for \(y\); in GSL terms, \(a_{y} = 1\) and \(a_{dy/dt} = 0\); in ARKode terms, the \(\eta^{abs}_i\) are identical
+                 | X'    Double Double -- ^ absolute and relative tolerance for \(\dot{y}\); in GSL terms, \(a_{y} = 0\) and \(a_{dy/dt} = 1\); in ARKode terms, the latter is treated as the relative tolerance for \(y\) so this is the same as specifying 'X' which may be entirely incorrect for the given problem
+                 | XX'   Double Double Double Double -- ^ include both via relative tolerance
+                                                     -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
+                 | ScXX' Double Double Double Double (Vector Double) -- ^ scale absolute tolerance of \(y_i\); in ARKode terms, \(a_{{dy}/{dt}}\) is ignored, \(\eta^{abs}_i = s_i \epsilon^{abs}\) and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
diff --git a/packages/sundials/src/Numeric/Sundials/ODEOpts.hs b/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
new file mode 100644
index 0000000..027d99a
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/ODEOpts.hs
@@ -0,0 +1,32 @@
+module Numeric.Sundials.ODEOpts where
+
+import           Data.Word (Word32)
+import qualified Data.Vector.Storable as VS
+
+import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix)
+
+
+type Jacobian = Double -> Vector Double -> Matrix Double
+
+data ODEOpts = ODEOpts {
+    maxNumSteps :: Word32
+  , minStep     :: Double
+  , relTol      :: Double
+  , absTols     :: VS.Vector Double
+  , initStep    :: Maybe Double
+  , maxFail     :: Word32
+  } deriving (Read, Show, Eq, Ord)
+
+data SundialsDiagnostics = SundialsDiagnostics {
+    aRKodeGetNumSteps               :: Int
+  , aRKodeGetNumStepAttempts        :: Int
+  , aRKodeGetNumRhsEvals_fe         :: Int
+  , aRKodeGetNumRhsEvals_fi         :: Int
+  , aRKodeGetNumLinSolvSetups       :: Int
+  , aRKodeGetNumErrTestFails        :: Int
+  , aRKodeGetNumNonlinSolvIters     :: Int
+  , aRKodeGetNumNonlinSolvConvFails :: Int
+  , aRKDlsGetNumJacEvals            :: Int
+  , aRKDlsGetNumRhsEvals            :: Int
+  } deriving Show
+
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c
new file mode 100644
index 0000000..f0ca592
--- /dev/null
+++ b/packages/sundials/src/helpers.c
@@ -0,0 +1,44 @@
+#include <stdio.h>
+#include <math.h>
+#include <arkode/arkode.h>                 /* prototypes for ARKODE fcts., consts. */
+#include <nvector/nvector_serial.h>        /* serial N_Vector types, fcts., macros */
+#include <sunmatrix/sunmatrix_dense.h>     /* access to dense SUNMatrix            */
+#include <sunlinsol/sunlinsol_dense.h>     /* access to dense SUNLinearSolver      */
+#include <arkode/arkode_direct.h>          /* access to ARKDls interface           */
+#include <sundials/sundials_types.h>       /* definition of type realtype          */
+#include <sundials/sundials_math.h>
+
+/* Check function return value...
+    opt == 0 means SUNDIALS function allocates memory so check if
+             returned NULL pointer
+    opt == 1 means SUNDIALS function returns a flag so check if
+             flag >= 0
+    opt == 2 means function allocates memory so check if returned
+             NULL pointer  
+*/
+int check_flag(void *flagvalue, const char *funcname, int opt)
+{
+  int *errflag;
+
+  /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
+  if (opt == 0 && flagvalue == NULL) {
+    fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
+            funcname);
+    return 1; }
+
+  /* Check if flag < 0 */
+  else if (opt == 1) {
+    errflag = (int *) flagvalue;
+    if (*errflag < 0) {
+      fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
+              funcname, *errflag);
+      return 1; }}
+
+  /* Check if function returned NULL pointer - no memory allocated */
+  else if (opt == 2 && flagvalue == NULL) {
+    fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
+            funcname);
+    return 1; }
+
+  return 0;
+}
diff --git a/packages/sundials/src/helpers.h b/packages/sundials/src/helpers.h
new file mode 100644
index 0000000..3d8fbc0
--- /dev/null
+++ b/packages/sundials/src/helpers.h
@@ -0,0 +1,9 @@
+/* Check function return value...
+    opt == 0 means SUNDIALS function allocates memory so check if
+             returned NULL pointer
+    opt == 1 means SUNDIALS function returns a flag so check if
+             flag >= 0
+    opt == 2 means function allocates memory so check if returned
+             NULL pointer
+*/
+int check_flag(void *flagvalue, const char *funcname, int opt);
diff --git a/packages/tests/hmatrix-tests.cabal b/packages/tests/hmatrix-tests.cabal
index 00f3a38..31fa32e 100644
--- a/packages/tests/hmatrix-tests.cabal
+++ b/packages/tests/hmatrix-tests.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-tests
-Version:            0.6.0.0
+Version:            0.19.0.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
diff --git a/packages/tests/src/Numeric/GSL/Tests.hs b/packages/tests/src/Numeric/GSL/Tests.hs
index 025427b..ed15935 100644
--- a/packages/tests/src/Numeric/GSL/Tests.hs
+++ b/packages/tests/src/Numeric/GSL/Tests.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-}
+{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns -fno-warn-missing-signatures #-}
 {- |
 Module      :  Numeric.GLS.Tests
 Copyright   :  (c) Alberto Ruiz 2014
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 2aefc87..4ed1462 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -1,5 +1,5 @@
 {-# LANGUAGE CPP #-}
-{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-}
+{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns -fno-warn-missing-signatures #-}
 {-# LANGUAGE DataKinds #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE FlexibleContexts #-}
@@ -31,7 +31,7 @@ module Numeric.LinearAlgebra.Tests(
 --, runBigTests
 ) where
 
-import Numeric.LinearAlgebra hiding (unitary)
+import Numeric.LinearAlgebra
 import Numeric.LinearAlgebra.Devel
 import Numeric.LinearAlgebra.Static(L)
 import Numeric.LinearAlgebra.Tests.Instances
@@ -514,7 +514,7 @@ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
 
 --------------------------------------------------------------------------------
 
-sliceTest = utest "slice test" $ and
+_sliceTest = TestList
     [ testSlice (chol . trustSym)  (gen 5 :: Matrix R)
     , testSlice (chol . trustSym)  (gen 5 :: Matrix C)
     , testSlice qr    (rec :: Matrix R)
@@ -617,7 +617,7 @@ sliceTest = utest "slice test" $ and
 
     test_qrgr n t x = qrgr n (QR x t)
 
-    ok_qrgr x = simeq 1E-15 q q'
+    ok_qrgr x = TestCase . assertBool "ok_qrgr" $ simeq 1E-15 q q'
       where
         (q,_) = qr x
         atau = qrRaw x
@@ -646,7 +646,8 @@ sliceTest = utest "slice test" $ and
     rec :: Numeric t => Matrix t
     rec = subMatrix (0,0) (4,5) (gen 5)
 
-    testSlice f x@(size->sz@(r,c)) = all (==f x) (map f (g y1 ++ g y2))
+    testSlice f x@(size->sz@(r,c)) =
+      TestList . map (TestCase . assertEqual "" (f x)) $ (map f (g y1 ++ g y2))
       where
         subm = subMatrix
         g y = [ subm (a*r,b*c) sz y | a <-[0..2], b <- [0..2]]
@@ -841,7 +842,7 @@ runTests n = do
         , staticTest
         , intTest
         , modularTest
-        , sliceTest
+        -- , sliceTest
         ]
     when (errors c + failures c > 0) exitFailure
     return ()
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index f0bddd0..59230e0 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,4 +1,8 @@
 {-# LANGUAGE CPP, FlexibleContexts, UndecidableInstances, FlexibleInstances, ScopedTypeVariables #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Tests.Instances
@@ -62,7 +66,7 @@ instance KnownNat n => Arbitrary (Static.R n) where
         n :: Int
         n = fromIntegral (natVal (Proxy :: Proxy n))
 
-    shrink v = []
+    shrink _v = []
 
 instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
     arbitrary = do
@@ -89,7 +93,7 @@ instance (KnownNat n, KnownNat m) => Arbitrary (Static.L m n) where
         n :: Int
         n = fromIntegral (natVal (Proxy :: Proxy n))
 
-    shrink mat = []
+    shrink _mat = []
 
 -- a square matrix
 newtype (Sq a) = Sq (Matrix a) deriving Show
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index e3a6242..6cd3a9e 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -3,6 +3,8 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE DataKinds #-}
 
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Tests.Properties
@@ -51,14 +53,13 @@ module Numeric.LinearAlgebra.Tests.Properties (
   , staticVectorBinaryFailProp
 ) where
 
-import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)
+import Numeric.LinearAlgebra.HMatrix hiding (Testable)
 import qualified Numeric.LinearAlgebra.Static as Static
 import Test.QuickCheck
 
 import Data.Binary
 import Data.Binary.Get (runGet)
 import Data.Either (isLeft)
-import Debug.Trace (traceShowId)
 #if MIN_VERSION_base(4,11,0)
 import Prelude hiding ((<>))
 #endif
diff --git a/packages/tests/src/TestBase.hs b/packages/tests/src/TestBase.hs
index 23fd675..51867b1 100644
--- a/packages/tests/src/TestBase.hs
+++ b/packages/tests/src/TestBase.hs
@@ -1,3 +1,4 @@
 import Numeric.LinearAlgebra.Tests
 
+main :: IO ()
 main = runTests 20
diff --git a/packages/tests/src/TestGSL.hs b/packages/tests/src/TestGSL.hs
index 112422d..cc6b1e7 100644
--- a/packages/tests/src/TestGSL.hs
+++ b/packages/tests/src/TestGSL.hs
@@ -1,3 +1,4 @@
 import Numeric.GSL.Tests
 
+main :: IO ()
 main = runTests 20
diff --git a/shell.nix b/shell.nix
new file mode 100644
index 0000000..42bbe08
--- /dev/null
+++ b/shell.nix
@@ -0,0 +1,26 @@
+{ nixpkgs ? import <nixpkgs> {} }:
+let p = nixpkgs; in
+
+  let
+  myStack = p.writeShellScriptBin "stack" ''
+    exec ${p.stack}/bin/stack --no-nix --no-docker --system-ghc $STACK_IN_NIX_EXTRA_ARGS "$@"
+  '';
+in
+p.haskell.lib.buildStackProject {
+  name = "hmatrix";
+
+  buildInputs = [
+    p.zlib
+    p.sundials
+    p.blas
+    p.gfortran.cc
+    p.liblapack
+    p.gsl
+    p.glpk
+    p.pkgconfig
+    p.stack
+    p.sundials
+    myStack
+  ] ++ (if p.stdenv.isDarwin then [p.darwin.apple_sdk.frameworks.Accelerate] else [])
+    ++ (if p.stdenv.isDarwin then [p.darwin.apple_sdk.frameworks.Cocoa] else []);
+}
diff --git a/stack.yaml b/stack.yaml
index 1198e57..9b858c0 100644
--- a/stack.yaml
+++ b/stack.yaml
@@ -13,5 +13,12 @@ packages:
 - packages/gsl/
 - packages/glpk/
 - packages/base/
-extra-deps: []
-resolver: lts-8.15
+- packages/sundials/
+- examples/
+extra-deps:
+- diagrams-rasterific-1.4
+- plots-0.1.0.2
+resolver: lts-10.9
+nix:
+  path: [nixpkgs=./nixpkgs.nix]
+  shell-file: shell.nix