67 files changed, 1664 insertions, 62 deletions
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 d64523d..46e9cf3 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 ComplexOf x = Complex (RealOf 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 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
+Version:            0.19.0.0
-License:            GPL
+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
+Version:            0.19.0.0
-License:            GPL
-License-file:       LICENSE
-Author:             Alberto Ruiz
-Maintainer:         Alberto Ruiz <aruiz@um.es>
-Stability:          provisional
-Homepage:           https://github.com/albertoruiz/hmatrix
 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
+    let v g = g (fi $ V.length x) p
-        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
+Version:            0.19.0.0
-License:            GPL
+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
Binary files differ
diff --git a/packages/sundials/hmatrix-sundials.cabal b/packages/sundials/hmatrix-sundials.cabal
new file mode 100644
index 0000000..388f1db
--- /dev/null
+++ b/packages/sundials/hmatrix-sundials.cabal
@@ -0,0 +1,57 @@
+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
+  other-extensions:    QuasiQuotes
+  hs-source-dirs:      src
+  exposed-modules:     Numeric.Sundials.ARKode.ODE
+  other-modules:       Types,
+                       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:       Types,
+                       Numeric.Sundials.ARKode.ODE,
+                       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
+  c-sources:           src/helpers.c src/helpers.h
+  default-language:    Haskell2010
diff --git a/packages/sundials/src/Arkode.hsc b/packages/sundials/src/Arkode.hsc
new file mode 100644
index 0000000..9db37b5
--- /dev/null
+++ b/packages/sundials/src/Arkode.hsc
@@ -0,0 +1,114 @@
+module Arkode where
+import Foreign
+import Foreign.C.Types
+#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>
+#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
+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/Main.hs b/packages/sundials/src/Main.hs
new file mode 100644
index 0000000..729d35a
--- /dev/null
+++ b/packages/sundials/src/Main.hs
@@ -0,0 +1,138 @@
+{-# OPTIONS_GHC -Wall #-}
+import           Numeric.Sundials.ARKode.ODE
+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
+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 = 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 = 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 = 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 = odeSolveV (SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
+  let res2b = odeSolveV (TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
+  let maxDiff = maximum $ map abs $
+                zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)
+  hspec $ describe "Compare results" $ do
+    it "for two different RK methods" $
+      maxDiff < 1.0e-6
+  let res3 = 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))
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..e5a2e4d
--- /dev/null
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -0,0 +1,898 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# 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.
+--
+-- <https://computation.llnl.gov/projects/sundials/sundials-software>
+--
+-- 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))
+-- @
+--
+-- 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(..)
+                                   , Jacobian
+                                   , SundialsDiagnostics(..)
+                                   ) 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
+import           Foreign.Ptr (Ptr)
+import           Foreign.ForeignPtr (newForeignPtr_)
+import           Foreign.Storable (Storable)
+import qualified Data.Vector.Storable as V
+import qualified Data.Vector.Storable.Mutable as VM
+import           Data.Coerce (coerce)
+import           System.IO.Unsafe (unsafePerformIO)
+import           GHC.Generics
+import           Numeric.LinearAlgebra.Devel (createVector)
+import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
+                                                subMatrix, rows, cols, toLists,
+                                                size, subVector)
+import qualified Types as T
+import           Arkode
+import qualified Arkode as B
+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 "Arkode_hsc.h"
+getDataFromContents :: Int -> Ptr T.SunVector -> IO (V.Vector CDouble)
+getDataFromContents len ptr = do
+  qtr <- B.getContentPtr ptr
+  rtr <- B.getData qtr
+  vectorFromC len rtr
+-- FIXME: Potentially an instance of Storable
+_getMatrixDataFromContents :: Ptr T.SunMatrix -> IO T.SunMatrix
+_getMatrixDataFromContents ptr = do
+  qtr <- B.getContentMatrixPtr ptr
+  rs  <- B.getNRows qtr
+  cs  <- B.getNCols qtr
+  rtr <- B.getMatrixData qtr
+  vs  <- vectorFromC (fromIntegral $ rs * cs) rtr
+  return $ T.SunMatrix { T.rows = rs, T.cols = cs, T.vals = vs }
+putMatrixDataFromContents :: T.SunMatrix -> Ptr T.SunMatrix -> IO ()
+putMatrixDataFromContents mat ptr = do
+  let rs = T.rows mat
+      cs = T.cols mat
+      vs = T.vals mat
+  qtr <- B.getContentMatrixPtr ptr
+  B.putNRows rs qtr
+  B.putNCols cs qtr
+  rtr <- B.getMatrixData qtr
+  vectorToC vs (fromIntegral $ rs * cs) rtr
+-- FIXME: END
+putDataInContents :: Storable a => V.Vector a -> Int -> Ptr b -> IO ()
+putDataInContents vec len ptr = do
+  qtr <- B.getContentPtr ptr
+  rtr <- B.getData qtr
+  vectorToC vec len rtr
+-- Utils
+vectorFromC :: Storable a => Int -> Ptr a -> IO (V.Vector a)
+vectorFromC len ptr = do
+  ptr' <- newForeignPtr_ ptr
+  V.freeze $ VM.unsafeFromForeignPtr0 ptr' len
+vectorToC :: Storable a => V.Vector a -> Int -> Ptr a -> IO ()
+vectorToC vec len ptr = do
+  ptr' <- newForeignPtr_ ptr
+  V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
+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
+type Jacobian = Double -> Vector Double -> Matrix Double
+-- | 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 =
+  case odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts of
+    Left c        -> error $ show c -- FIXME
+    -- FIXME: Can we do better than using lists?
+    Right (v, _d) -> (nR >< nC) (V.toList v)
+  where
+    us = toList ts
+    nR = length us
+    nC = size y0
+    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
+  case odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1)  Nothing g (V.fromList y0) (V.fromList $ toList ts) of
+    Left c        -> error $ show c -- FIXME
+    Right (v, _d) -> (nR >< nC) (V.toList v)
+  where
+    us = toList ts
+    nR = length us
+    nC = length y0
+    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 method control initStepSize f y0 tt of
+    Left c        -> error $ show c -- FIXME
+    Right (v, _d) -> (nR >< nC) (V.toList v)
+  where
+    nR = V.length tt
+    nC = V.length y0
+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
+  -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution
+odeSolveVWith method control initStepSize f y0 tt =
+  case solveOdeC (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
+                 (coerce f) (coerce y0) (coerce tt) of
+    Left c -> Left $ fromIntegral c
+    Right (v, d) -> Right (coerce v, d)
+  where
+    l = size y0
+    scise (X absTol relTol)                          = coerce (V.replicate l absTol, relTol)
+    scise (X' absTol relTol)                         = coerce (V.replicate l absTol, relTol)
+    scise (XX' absTol relTol yScale _yDotScale)      = coerce (V.replicate l absTol, yScale * relTol)
+    -- FIXME; Should we check that the length of ss is correct?
+    scise (ScXX' absTol relTol yScale _yDotScale ss) = coerce (V.map (* absTol) ss, yScale * relTol)
+    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 ->
+  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 method initStepSize jacH (absTols, relTol) 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: fMut is not actually mutatated
+  fMut <- V.thaw f0
+  tMut <- V.thaw 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 'Types'. 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
+                          putMatrixDataFromContents j jacS
+                          -- 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 *tMut))[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 *fMut))[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 *absTols))[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 Types                             */
+                         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;
+                         }
+                         /* FIXME: A hack for initial testing */
+                         flag = ARKodeSetMinStep(arkode_mem, 1.0e-12);
+                         if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1;
+                         flag = ARKodeSetMaxNumSteps(arkode_mem, 10000);
+                         if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1;
+                         /* Set routines */
+                         flag = ARKodeSVtolerances(arkode_mem, $(double relTol), 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 *tMut))[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) $ (B.arkSMax >< B.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 (B.arkSMax * B.arkSMax)
+  btAsMut <- V.thaw btAs
+  btCs  :: V.Vector CDouble <- createVector B.arkSMax
+  btBs  :: V.Vector CDouble <- createVector B.arkSMax
+  btB2s :: V.Vector CDouble <- createVector B.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/Types.hs b/packages/sundials/src/Types.hs
new file mode 100644
index 0000000..04e4280
--- /dev/null
+++ b/packages/sundials/src/Types.hs
@@ -0,0 +1,40 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE MultiWayIf #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE EmptyDataDecls #-}
+module Types where
+import           Foreign.C.Types
+import qualified Language.Haskell.TH as TH
+import qualified Language.C.Types as CT
+import qualified Data.Map as Map
+import           Language.C.Inline.Context
+import qualified Data.Vector.Storable as V
+data SunVector
+data SunMatrix = SunMatrix { rows :: CInt
+                           , cols :: CInt
+                           , vals :: V.Vector CDouble
+                           }
+-- FIXME: Is this true?
+type SunIndexType = CLong
+sunTypesTable :: Map.Map CT.TypeSpecifier TH.TypeQ
+sunTypesTable = Map.fromList
+  [
+    (CT.TypeName "sunindextype", [t| SunIndexType |] )
+  , (CT.TypeName "SunVector", [t| SunVector |] )
+  , (CT.TypeName "SunMatrix", [t| SunMatrix |] )
+  ]
+sunCtx :: Context
+sunCtx = mempty {ctxTypesTable = sunTypesTable}
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