summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
Diffstat (limited to 'packages')
-rw-r--r--packages/base/CHANGELOG7
-rw-r--r--packages/base/hmatrix.cabal2
-rw-r--r--packages/base/src/Internal/Algorithms.hs2
-rw-r--r--packages/base/src/Internal/CG.hs2
-rw-r--r--packages/base/src/Internal/Chain.hs2
-rw-r--r--packages/base/src/Internal/Container.hs2
-rw-r--r--packages/base/src/Internal/Devel.hs1
-rw-r--r--packages/base/src/Internal/Element.hs23
-rw-r--r--packages/base/src/Internal/IO.hs16
-rw-r--r--packages/base/src/Internal/LAPACK.hs2
-rw-r--r--packages/base/src/Internal/Matrix.hs87
-rw-r--r--packages/base/src/Internal/Modular.hs3
-rw-r--r--packages/base/src/Internal/Numeric.hs10
-rw-r--r--packages/base/src/Internal/ST.hs12
-rw-r--r--packages/base/src/Internal/Sparse.hs2
-rw-r--r--packages/base/src/Internal/Static.hs3
-rw-r--r--packages/base/src/Internal/Util.hs2
-rw-r--r--packages/base/src/Internal/Vector.hs7
-rw-r--r--packages/base/src/Internal/Vectorized.hs36
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/HMatrix.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Static.hs2
-rw-r--r--packages/base/src/Numeric/Matrix.hs14
-rw-r--r--packages/base/src/Numeric/Vector.hs14
-rw-r--r--packages/glpk/hmatrix-glpk.cabal4
-rw-r--r--packages/gsl/CHANGELOG3
-rw-r--r--packages/gsl/hmatrix-gsl.cabal26
-rw-r--r--packages/gsl/src/Graphics/Plot.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Fitting.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Fourier.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Integration.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Internal.hs6
-rw-r--r--packages/gsl/src/Numeric/GSL/Interpolation.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/LinearAlgebra.hs3
-rw-r--r--packages/gsl/src/Numeric/GSL/Minimization.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/ODE.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Root.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Vector.hs3
-rw-r--r--packages/sparse/hmatrix-sparse.cabal4
-rw-r--r--packages/special/hmatrix-special.cabal4
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Bessel.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Coulomb.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Coupling.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Exp.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Gamma.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Gegenbauer.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Legendre.hs3
-rw-r--r--packages/special/lib/Numeric/GSL/Special/Trig.hs4
-rw-r--r--packages/sundials/ChangeLog.md5
-rw-r--r--packages/sundials/LICENSE30
-rw-r--r--packages/sundials/README.md8
-rw-r--r--packages/sundials/Setup.hs2
-rw-r--r--packages/sundials/diagrams/brusselator.pngbin0 -> 27362 bytes
-rw-r--r--packages/sundials/hmatrix-sundials.cabal57
-rw-r--r--packages/sundials/src/Arkode.hsc114
-rw-r--r--packages/sundials/src/Main.hs138
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs898
-rw-r--r--packages/sundials/src/Types.hs40
-rw-r--r--packages/sundials/src/helpers.c44
-rw-r--r--packages/sundials/src/helpers.h9
-rw-r--r--packages/tests/hmatrix-tests.cabal2
-rw-r--r--packages/tests/src/Numeric/GSL/Tests.hs2
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs13
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs8
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs5
-rw-r--r--packages/tests/src/TestBase.hs1
-rw-r--r--packages/tests/src/TestGSL.hs1
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 @@
2-------- 2--------
3 3
4 * Many new functions and instances in the Static module 4 * Many new functions and instances in the Static module
5 5
6 * meanCov and gaussianSample use Herm type 6 * meanCov and gaussianSample use Herm type
7 7
8 * thinQR, thinRQ 8 * thinQR, thinRQ
9 9
10 * compactSVDTol 10 * compactSVDTol
11 11
12 * unitary changed to normalize, also admits Vector (Complex Double) 12 * unitary changed to normalize, also admits Vector (Complex Double)
13 13
140.17.0.0 140.17.0.0
15-------- 15--------
16 16
@@ -288,4 +288,3 @@
288 * added NFData instances for Matrix and Vector. 288 * added NFData instances for Matrix and Vector.
289 289
290 * liftVector, liftVector2 replaced by mapVector, zipVector. 290 * liftVector, liftVector2 replaced by mapVector, zipVector.
291
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 @@
1Name: hmatrix 1Name: hmatrix
2Version: 0.18.2.0 2Version: 0.19.0.0
3License: BSD3 3License: BSD3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: 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 @@
4{-# LANGUAGE UndecidableInstances #-} 4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE TypeFamilies #-} 5{-# LANGUAGE TypeFamilies #-}
6 6
7{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
8
7----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
8{- | 10{- |
9Module : Internal.Algorithms 11Module : 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 @@
1{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} 1{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
2{-# LANGUAGE RecordWildCards #-} 2{-# LANGUAGE RecordWildCards #-}
3 3
4{-# OPTIONS_GHC -fno-warn-orphans #-}
5
4module Internal.CG( 6module Internal.CG(
5 cgSolve, cgSolve', 7 cgSolve, cgSolve',
6 CGState(..), R, V 8 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 @@
1{-# LANGUAGE FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2 2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4
3----------------------------------------------------------------------------- 5-----------------------------------------------------------------------------
4-- | 6-- |
5-- Module : Internal.Chain 7-- 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 @@
5{-# LANGUAGE FunctionalDependencies #-} 5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-} 6{-# LANGUAGE UndecidableInstances #-}
7 7
8{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}
9
8----------------------------------------------------------------------------- 10-----------------------------------------------------------------------------
9-- | 11-- |
10-- Module : Internal.Container 12-- 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
54 54
55-- | postfix error code check 55-- | postfix error code check
56infixl 0 #| 56infixl 0 #|
57(#|) :: IO CInt -> String -> IO ()
57(#|) = flip check 58(#|) = flip check
58 59
59-- | Error capture and conversion to Maybe 60-- | 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 @@
4{-# LANGUAGE UndecidableInstances #-} 4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-} 5{-# LANGUAGE MultiParamTypeClasses #-}
6 6
7{-# OPTIONS_GHC -fno-warn-orphans #-}
8
7----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
8-- | 10-- |
9-- Module : Data.Packed.Matrix 11-- Module : Data.Packed.Matrix
@@ -31,6 +33,7 @@ import Data.List.Split(chunksOf)
31import Foreign.Storable(Storable) 33import Foreign.Storable(Storable)
32import System.IO.Unsafe(unsafePerformIO) 34import System.IO.Unsafe(unsafePerformIO)
33import Control.Monad(liftM) 35import Control.Monad(liftM)
36import Foreign.C.Types(CInt)
34 37
35------------------------------------------------------------------- 38-------------------------------------------------------------------
36 39
@@ -53,8 +56,10 @@ instance (Show a, Element a) => (Show (Matrix a)) where
53 show m | rows m == 0 || cols m == 0 = sizes m ++" []" 56 show m | rows m == 0 || cols m == 0 = sizes m ++" []"
54 show m = (sizes m++) . dsp . map (map show) . toLists $ m 57 show m = (sizes m++) . dsp . map (map show) . toLists $ m
55 58
59sizes :: Matrix t -> [Char]
56sizes m = "("++show (rows m)++"><"++show (cols m)++")\n" 60sizes m = "("++show (rows m)++"><"++show (cols m)++")\n"
57 61
62dsp :: [[[Char]]] -> [Char]
58dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp 63dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
59 where 64 where
60 mt = transpose as 65 mt = transpose as
@@ -73,6 +78,7 @@ instance (Element a, Read a) => Read (Matrix a) where
73 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims 78 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
74 79
75 80
81breakAt :: Eq a => a -> [a] -> ([a], [a])
76breakAt c l = (a++[c],tail b) where 82breakAt c l = (a++[c],tail b) where
77 (a,b) = break (==c) l 83 (a,b) = break (==c) l
78 84
@@ -88,7 +94,8 @@ data Extractor
88 | Drop Int 94 | Drop Int
89 | DropLast Int 95 | DropLast Int
90 deriving Show 96 deriving Show
91 97
98ppext :: Extractor -> [Char]
92ppext All = ":" 99ppext All = ":"
93ppext (Range a 1 c) = printf "%d:%d" a c 100ppext (Range a 1 c) = printf "%d:%d" a c
94ppext (Range a b c) = printf "%d:%d:%d" a b c 101ppext (Range a b c) = printf "%d:%d:%d" a b c
@@ -128,10 +135,14 @@ ppext (DropLast n) = printf "DropLast %d" n
128infixl 9 ?? 135infixl 9 ??
129(??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t 136(??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t
130 137
138minEl :: Vector CInt -> CInt
131minEl = toScalarI Min 139minEl = toScalarI Min
140maxEl :: Vector CInt -> CInt
132maxEl = toScalarI Max 141maxEl = toScalarI Max
142cmodi :: Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt -> Vector Foreign.C.Types.CInt
133cmodi = vectorMapValI ModVS 143cmodi = vectorMapValI ModVS
134 144
145extractError :: Matrix t1 -> (Extractor, Extractor) -> t
135extractError m (e1,e2)= error $ printf "can't extract (%s,%s) from matrix %dx%d" (ppext e1::String) (ppext e2::String) (rows m) (cols m) 146extractError m (e1,e2)= error $ printf "can't extract (%s,%s) from matrix %dx%d" (ppext e1::String) (ppext e2::String) (rows m) (cols m)
136 147
137m ?? (Range a s b,e) | s /= 1 = m ?? (Pos (idxs [a,a+s .. b]), e) 148m ?? (Range a s b,e) | s /= 1 = m ?? (Pos (idxs [a,a+s .. b]), e)
@@ -232,8 +243,10 @@ disp = putStr . dispf 2
232fromBlocks :: Element t => [[Matrix t]] -> Matrix t 243fromBlocks :: Element t => [[Matrix t]] -> Matrix t
233fromBlocks = fromBlocksRaw . adaptBlocks 244fromBlocks = fromBlocksRaw . adaptBlocks
234 245
246fromBlocksRaw :: Element t => [[Matrix t]] -> Matrix t
235fromBlocksRaw mms = joinVert . map joinHoriz $ mms 247fromBlocksRaw mms = joinVert . map joinHoriz $ mms
236 248
249adaptBlocks :: Element t => [[Matrix t]] -> [[Matrix t]]
237adaptBlocks ms = ms' where 250adaptBlocks ms = ms' where
238 bc = case common length ms of 251 bc = case common length ms of
239 Just c -> c 252 Just c -> c
@@ -486,6 +499,9 @@ liftMatrix2Auto f m1 m2
486 m2' = conformMTo (r,c) m2 499 m2' = conformMTo (r,c) m2
487 500
488-- FIXME do not flatten if equal order 501-- FIXME do not flatten if equal order
502lM :: (Storable t, Element t1, Element t2)
503 => (Vector t1 -> Vector t2 -> Vector t)
504 -> Matrix t1 -> Matrix t2 -> Matrix t
489lM f m1 m2 = matrixFromVector 505lM f m1 m2 = matrixFromVector
490 RowMajor 506 RowMajor
491 (max' (rows m1) (rows m2)) 507 (max' (rows m1) (rows m2))
@@ -504,6 +520,7 @@ compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
504 520
505------------------------------------------------------------ 521------------------------------------------------------------
506 522
523toBlockRows :: Element t => [Int] -> Matrix t -> [Matrix t]
507toBlockRows [r] m 524toBlockRows [r] m
508 | r == rows m = [m] 525 | r == rows m = [m]
509toBlockRows rs m 526toBlockRows rs m
@@ -513,6 +530,7 @@ toBlockRows rs m
513 szs = map (* cols m) rs 530 szs = map (* cols m) rs
514 g k = (k><0)[] 531 g k = (k><0)[]
515 532
533toBlockCols :: Element t => [Int] -> Matrix t -> [Matrix t]
516toBlockCols [c] m | c == cols m = [m] 534toBlockCols [c] m | c == cols m = [m]
517toBlockCols cs m = map trans . toBlockRows cs . trans $ m 535toBlockCols cs m = map trans . toBlockRows cs . trans $ m
518 536
@@ -576,7 +594,7 @@ Just (3><3)
576mapMatrixWithIndexM 594mapMatrixWithIndexM
577 :: (Element a, Storable b, Monad m) => 595 :: (Element a, Storable b, Monad m) =>
578 ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b) 596 ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
579mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m 597mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m
580 where 598 where
581 c = cols m 599 c = cols m
582 600
@@ -598,4 +616,3 @@ mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
598 616
599mapMatrix :: (Element a, Element b) => (a -> b) -> Matrix a -> Matrix b 617mapMatrix :: (Element a, Element b) => (a -> b) -> Matrix a -> Matrix b
600mapMatrix f = liftMatrix (mapVector f) 618mapMatrix f = liftMatrix (mapVector f)
601
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
20import Internal.Vector 20import Internal.Vector
21import Internal.Matrix 21import Internal.Matrix
22import Internal.Vectorized 22import Internal.Vectorized
23import Text.Printf(printf) 23import Text.Printf(printf, PrintfArg, PrintfType)
24import Data.List(intersperse,transpose) 24import Data.List(intersperse,transpose)
25import Data.Complex 25import Data.Complex
26 26
@@ -78,12 +78,18 @@ disps d x = sdims x ++ " " ++ formatScaled d x
78dispf :: Int -> Matrix Double -> String 78dispf :: Int -> Matrix Double -> String
79dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x 79dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
80 80
81sdims :: Matrix t -> [Char]
81sdims x = show (rows x) ++ "x" ++ show (cols x) 82sdims x = show (rows x) ++ "x" ++ show (cols x)
82 83
84formatFixed :: (Show a, Text.Printf.PrintfArg t, Element t)
85 => a -> Matrix t -> String
83formatFixed d x = format " " (printf ("%."++show d++"f")) $ x 86formatFixed d x = format " " (printf ("%."++show d++"f")) $ x
84 87
88isInt :: Matrix Double -> Bool
85isInt = all lookslikeInt . toList . flatten 89isInt = all lookslikeInt . toList . flatten
86 90
91formatScaled :: (Text.Printf.PrintfArg b, RealFrac b, Floating b, Num t, Element b, Show t)
92 => t -> Matrix b -> [Char]
87formatScaled dec t = "E"++show o++"\n" ++ ss 93formatScaled dec t = "E"++show o++"\n" ++ ss
88 where ss = format " " (printf fmt. g) t 94 where ss = format " " (printf fmt. g) t
89 g x | o >= 0 = x/10^(o::Int) 95 g x | o >= 0 = x/10^(o::Int)
@@ -133,14 +139,18 @@ showComplex d (a:+b)
133 s2 = if b<0 then "-" else "" 139 s2 = if b<0 then "-" else ""
134 s3 = if b<0 then "-" else "+" 140 s3 = if b<0 then "-" else "+"
135 141
142shcr :: (Show a, Show t1, Text.Printf.PrintfType t, Text.Printf.PrintfArg t1, RealFrac t1)
143 => a -> t1 -> t
136shcr d a | lookslikeInt a = printf "%.0f" a 144shcr d a | lookslikeInt a = printf "%.0f" a
137 | otherwise = printf ("%."++show d++"f") a 145 | otherwise = printf ("%."++show d++"f") a
138 146
139 147lookslikeInt :: (Show a, RealFrac a) => a -> Bool
140lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx 148lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
141 where shx = show x 149 where shx = show x
142 150
151isZero :: Show a => a -> Bool
143isZero x = show x `elem` ["0.0","-0.0"] 152isZero x = show x `elem` ["0.0","-0.0"]
153isOne :: Show a => a -> Bool
144isOne x = show x `elem` ["1.0","-1.0"] 154isOne x = show x `elem` ["1.0","-1.0"]
145 155
146-- | Pretty print a complex matrix with at most n decimal digits. 156-- | Pretty print a complex matrix with at most n decimal digits.
@@ -168,6 +178,6 @@ loadMatrix f = do
168 else 178 else
169 return (reshape c v) 179 return (reshape c v)
170 180
171 181loadMatrix' :: FilePath -> IO (Maybe (Matrix Double))
172loadMatrix' name = mbCatch (loadMatrix name) 182loadMatrix' name = mbCatch (loadMatrix name)
173 183
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 @@
1{-# LANGUAGE TypeOperators #-} 1{-# LANGUAGE TypeOperators #-}
2{-# LANGUAGE ViewPatterns #-} 2{-# LANGUAGE ViewPatterns #-}
3 3
4{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
5
4----------------------------------------------------------------------------- 6-----------------------------------------------------------------------------
5-- | 7-- |
6-- Module : Numeric.LinearAlgebra.LAPACK 8-- 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
57cols = icols 57cols = icols
58{-# INLINE cols #-} 58{-# INLINE cols #-}
59 59
60size :: Matrix t -> (Int, Int)
60size m = (irows m, icols m) 61size m = (irows m, icols m)
61{-# INLINE size #-} 62{-# INLINE size #-}
62 63
64rowOrder :: Matrix t -> Bool
63rowOrder m = xCol m == 1 || cols m == 1 65rowOrder m = xCol m == 1 || cols m == 1
64{-# INLINE rowOrder #-} 66{-# INLINE rowOrder #-}
65 67
68colOrder :: Matrix t -> Bool
66colOrder m = xRow m == 1 || rows m == 1 69colOrder m = xRow m == 1 || rows m == 1
67{-# INLINE colOrder #-} 70{-# INLINE colOrder #-}
68 71
72is1d :: Matrix t -> Bool
69is1d (size->(r,c)) = r==1 || c==1 73is1d (size->(r,c)) = r==1 || c==1
70{-# INLINE is1d #-} 74{-# INLINE is1d #-}
71 75
72-- data is not contiguous 76-- data is not contiguous
77isSlice :: Storable t => Matrix t -> Bool
73isSlice m@(size->(r,c)) = r*c < dim (xdat m) 78isSlice m@(size->(r,c)) = r*c < dim (xdat m)
74{-# INLINE isSlice #-} 79{-# INLINE isSlice #-}
75 80
@@ -136,16 +141,20 @@ instance Storable t => TransArray (Matrix t)
136 {-# INLINE applyRaw #-} 141 {-# INLINE applyRaw #-}
137 142
138infixr 1 # 143infixr 1 #
144(#) :: TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
139a # b = apply a b 145a # b = apply a b
140{-# INLINE (#) #-} 146{-# INLINE (#) #-}
141 147
148(#!) :: (TransArray c, TransArray c1) => c1 -> c -> Trans c1 (Trans c (IO r)) -> IO r
142a #! b = a # b # id 149a #! b = a # b # id
143{-# INLINE (#!) #-} 150{-# INLINE (#!) #-}
144 151
145-------------------------------------------------------------------------------- 152--------------------------------------------------------------------------------
146 153
154copy :: Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
147copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1]) 155copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])
148 156
157extractAll :: Element t => MatrixOrder -> Matrix t -> Matrix t
149extractAll ord m = unsafePerformIO (copy ord m) 158extractAll ord m = unsafePerformIO (copy ord m)
150 159
151{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. 160{- | 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)
224{-# INLINE (@@>) #-} 233{-# INLINE (@@>) #-}
225 234
226-- Unsafe matrix access without range checking 235-- Unsafe matrix access without range checking
236atM' :: Storable t => Matrix t -> Int -> Int -> t
227atM' m i j = xdat m `at'` (i * (xRow m) + j * (xCol m)) 237atM' m i j = xdat m `at'` (i * (xRow m) + j * (xCol m))
228{-# INLINE atM' #-} 238{-# INLINE atM' #-}
229 239
230------------------------------------------------------------------ 240------------------------------------------------------------------
231 241
242matrixFromVector :: Storable t => MatrixOrder -> Int -> Int -> Vector t -> Matrix t
232matrixFromVector _ 1 _ v@(dim->d) = Matrix { irows = 1, icols = d, xdat = v, xRow = d, xCol = 1 } 243matrixFromVector _ 1 _ v@(dim->d) = Matrix { irows = 1, icols = d, xdat = v, xRow = d, xCol = 1 }
233matrixFromVector _ _ 1 v@(dim->d) = Matrix { irows = d, icols = 1, xdat = v, xRow = 1, xCol = d } 244matrixFromVector _ _ 1 v@(dim->d) = Matrix { irows = d, icols = 1, xdat = v, xRow = 1, xCol = d }
234matrixFromVector o r c v 245matrixFromVector o r c v
@@ -388,18 +399,21 @@ subMatrix (r0,c0) (rt,ct) m
388 399
389-------------------------------------------------------------------------- 400--------------------------------------------------------------------------
390 401
402maxZ :: (Num t1, Ord t1, Foldable t) => t t1 -> t1
391maxZ xs = if minimum xs == 0 then 0 else maximum xs 403maxZ xs = if minimum xs == 0 then 0 else maximum xs
392 404
405conformMs :: Element t => [Matrix t] -> [Matrix t]
393conformMs ms = map (conformMTo (r,c)) ms 406conformMs ms = map (conformMTo (r,c)) ms
394 where 407 where
395 r = maxZ (map rows ms) 408 r = maxZ (map rows ms)
396 c = maxZ (map cols ms) 409 c = maxZ (map cols ms)
397 410
398 411conformVs :: Element t => [Vector t] -> [Vector t]
399conformVs vs = map (conformVTo n) vs 412conformVs vs = map (conformVTo n) vs
400 where 413 where
401 n = maxZ (map dim vs) 414 n = maxZ (map dim vs)
402 415
416conformMTo :: Element t => (Int, Int) -> Matrix t -> Matrix t
403conformMTo (r,c) m 417conformMTo (r,c) m
404 | size m == (r,c) = m 418 | size m == (r,c) = m
405 | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c)) 419 | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
@@ -407,18 +421,24 @@ conformMTo (r,c) m
407 | size m == (1,c) = repRows r m 421 | size m == (1,c) = repRows r m
408 | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c) 422 | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c)
409 423
424conformVTo :: Element t => Int -> Vector t -> Vector t
410conformVTo n v 425conformVTo n v
411 | dim v == n = v 426 | dim v == n = v
412 | dim v == 1 = constantD (v@>0) n 427 | dim v == 1 = constantD (v@>0) n
413 | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n 428 | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
414 429
430repRows :: Element t => Int -> Matrix t -> Matrix t
415repRows n x = fromRows (replicate n (flatten x)) 431repRows n x = fromRows (replicate n (flatten x))
432repCols :: Element t => Int -> Matrix t -> Matrix t
416repCols n x = fromColumns (replicate n (flatten x)) 433repCols n x = fromColumns (replicate n (flatten x))
417 434
435shSize :: Matrix t -> [Char]
418shSize = shDim . size 436shSize = shDim . size
419 437
438shDim :: (Show a, Show a1) => (a1, a) -> [Char]
420shDim (r,c) = "(" ++ show r ++"x"++ show c ++")" 439shDim (r,c) = "(" ++ show r ++"x"++ show c ++")"
421 440
441emptyM :: Storable t => Int -> Int -> Matrix t
422emptyM r c = matrixFromVector RowMajor r c (fromList[]) 442emptyM r c = matrixFromVector RowMajor r c (fromList[])
423 443
424---------------------------------------------------------------------- 444----------------------------------------------------------------------
@@ -433,6 +453,11 @@ instance (Storable t, NFData t) => NFData (Matrix t)
433 453
434--------------------------------------------------------------- 454---------------------------------------------------------------
435 455
456extractAux :: (Eq t3, Eq t2, TransArray c, Storable a, Storable t1,
457 Storable t, Num t3, Num t2, Integral t1, Integral t)
458 => (t3 -> t2 -> CInt -> Ptr t1 -> CInt -> Ptr t
459 -> Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt))
460 -> MatrixOrder -> c -> t3 -> Vector t1 -> t2 -> Vector t -> IO (Matrix a)
436extractAux f ord m moder vr modec vc = do 461extractAux f ord m moder vr modec vc = do
437 let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr 462 let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
438 nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc 463 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
452 477
453--------------------------------------------------------------- 478---------------------------------------------------------------
454 479
480setRectAux :: (TransArray c1, TransArray c)
481 => (CInt -> CInt -> Trans c1 (Trans c (IO CInt)))
482 -> Int -> Int -> c1 -> c -> IO ()
455setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect" 483setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"
456 484
457type SetRect x = I -> I -> x ::> x::> Ok 485type SetRect x = I -> I -> x ::> x::> Ok
@@ -465,19 +493,29 @@ foreign import ccall unsafe "setRectL" c_setRectL :: SetRect Z
465 493
466-------------------------------------------------------------------------------- 494--------------------------------------------------------------------------------
467 495
496sortG :: (Storable t, Storable a)
497 => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
468sortG f v = unsafePerformIO $ do 498sortG f v = unsafePerformIO $ do
469 r <- createVector (dim v) 499 r <- createVector (dim v)
470 (v #! r) f #|"sortG" 500 (v #! r) f #|"sortG"
471 return r 501 return r
472 502
503sortIdxD :: Vector Double -> Vector CInt
473sortIdxD = sortG c_sort_indexD 504sortIdxD = sortG c_sort_indexD
505sortIdxF :: Vector Float -> Vector CInt
474sortIdxF = sortG c_sort_indexF 506sortIdxF = sortG c_sort_indexF
507sortIdxI :: Vector CInt -> Vector CInt
475sortIdxI = sortG c_sort_indexI 508sortIdxI = sortG c_sort_indexI
509sortIdxL :: Vector Z -> Vector I
476sortIdxL = sortG c_sort_indexL 510sortIdxL = sortG c_sort_indexL
477 511
512sortValD :: Vector Double -> Vector Double
478sortValD = sortG c_sort_valD 513sortValD = sortG c_sort_valD
514sortValF :: Vector Float -> Vector Float
479sortValF = sortG c_sort_valF 515sortValF = sortG c_sort_valF
516sortValI :: Vector CInt -> Vector CInt
480sortValI = sortG c_sort_valI 517sortValI = sortG c_sort_valI
518sortValL :: Vector Z -> Vector Z
481sortValL = sortG c_sort_valL 519sortValL = sortG c_sort_valL
482 520
483foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt)) 521foreign 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
492 530
493-------------------------------------------------------------------------------- 531--------------------------------------------------------------------------------
494 532
533compareG :: (TransArray c, Storable t, Storable a)
534 => Trans c (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt)
535 -> c -> Vector t -> Vector a
495compareG f u v = unsafePerformIO $ do 536compareG f u v = unsafePerformIO $ do
496 r <- createVector (dim v) 537 r <- createVector (dim v)
497 (u # v #! r) f #|"compareG" 538 (u # v #! r) f #|"compareG"
498 return r 539 return r
499 540
541compareD :: Vector Double -> Vector Double -> Vector CInt
500compareD = compareG c_compareD 542compareD = compareG c_compareD
543compareF :: Vector Float -> Vector Float -> Vector CInt
501compareF = compareG c_compareF 544compareF = compareG c_compareF
545compareI :: Vector CInt -> Vector CInt -> Vector CInt
502compareI = compareG c_compareI 546compareI = compareG c_compareI
547compareL :: Vector Z -> Vector Z -> Vector CInt
503compareL = compareG c_compareL 548compareL = compareG c_compareL
504 549
505foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt))) 550foreign 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
509 554
510-------------------------------------------------------------------------------- 555--------------------------------------------------------------------------------
511 556
557selectG :: (TransArray c, TransArray c1, TransArray c2, Storable t, Storable a)
558 => Trans c2 (Trans c1 (CInt -> Ptr t -> Trans c (CInt -> Ptr a -> IO CInt)))
559 -> c2 -> c1 -> Vector t -> c -> Vector a
512selectG f c u v w = unsafePerformIO $ do 560selectG f c u v w = unsafePerformIO $ do
513 r <- createVector (dim v) 561 r <- createVector (dim v)
514 (c # u # v # w #! r) f #|"selectG" 562 (c # u # v # w #! r) f #|"selectG"
515 return r 563 return r
516 564
565selectD :: Vector CInt -> Vector Double -> Vector Double -> Vector Double -> Vector Double
517selectD = selectG c_selectD 566selectD = selectG c_selectD
567selectF :: Vector CInt -> Vector Float -> Vector Float -> Vector Float -> Vector Float
518selectF = selectG c_selectF 568selectF = selectG c_selectF
569selectI :: Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt
519selectI = selectG c_selectI 570selectI = selectG c_selectI
571selectL :: Vector CInt -> Vector Z -> Vector Z -> Vector Z -> Vector Z
520selectL = selectG c_selectL 572selectL = selectG c_selectL
573selectC :: Vector CInt
574 -> Vector (Complex Double)
575 -> Vector (Complex Double)
576 -> Vector (Complex Double)
577 -> Vector (Complex Double)
521selectC = selectG c_selectC 578selectC = selectG c_selectC
579selectQ :: Vector CInt
580 -> Vector (Complex Float)
581 -> Vector (Complex Float)
582 -> Vector (Complex Float)
583 -> Vector (Complex Float)
522selectQ = selectG c_selectQ 584selectQ = selectG c_selectQ
523 585
524type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt))))) 586type 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
532 594
533--------------------------------------------------------------------------- 595---------------------------------------------------------------------------
534 596
597remapG :: (TransArray c, TransArray c1, Storable t, Storable a)
598 => (CInt -> CInt -> CInt -> CInt -> Ptr t
599 -> Trans c1 (Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)))
600 -> Matrix t -> c1 -> c -> Matrix a
535remapG f i j m = unsafePerformIO $ do 601remapG f i j m = unsafePerformIO $ do
536 r <- createMatrix RowMajor (rows i) (cols i) 602 r <- createMatrix RowMajor (rows i) (cols i)
537 (i # j # m #! r) f #|"remapG" 603 (i # j # m #! r) f #|"remapG"
538 return r 604 return r
539 605
606remapD :: Matrix CInt -> Matrix CInt -> Matrix Double -> Matrix Double
540remapD = remapG c_remapD 607remapD = remapG c_remapD
608remapF :: Matrix CInt -> Matrix CInt -> Matrix Float -> Matrix Float
541remapF = remapG c_remapF 609remapF = remapG c_remapF
610remapI :: Matrix CInt -> Matrix CInt -> Matrix CInt -> Matrix CInt
542remapI = remapG c_remapI 611remapI = remapG c_remapI
612remapL :: Matrix CInt -> Matrix CInt -> Matrix Z -> Matrix Z
543remapL = remapG c_remapL 613remapL = remapG c_remapL
614remapC :: Matrix CInt
615 -> Matrix CInt
616 -> Matrix (Complex Double)
617 -> Matrix (Complex Double)
544remapC = remapG c_remapC 618remapC = remapG c_remapC
619remapQ :: Matrix CInt -> Matrix CInt -> Matrix (Complex Float) -> Matrix (Complex Float)
545remapQ = remapG c_remapQ 620remapQ = remapG c_remapQ
546 621
547type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt)))) 622type 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
555 630
556-------------------------------------------------------------------------------- 631--------------------------------------------------------------------------------
557 632
633rowOpAux :: (TransArray c, Storable a) =>
634 (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Trans c (IO CInt))
635 -> Int -> a -> Int -> Int -> Int -> Int -> c -> IO ()
558rowOpAux f c x i1 i2 j1 j2 m = do 636rowOpAux f c x i1 i2 j1 j2 m = do
559 px <- newArray [x] 637 px <- newArray [x]
560 (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp" 638 (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
573 651
574-------------------------------------------------------------------------------- 652--------------------------------------------------------------------------------
575 653
654gemmg :: (TransArray c1, TransArray c, TransArray c2, TransArray c3)
655 => Trans c3 (Trans c2 (Trans c1 (Trans c (IO CInt))))
656 -> c3 -> c2 -> c1 -> c -> IO ()
576gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg" 657gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"
577 658
578type Tgemm x = x :> x ::> x ::> x ::> Ok 659type Tgemm x = x :> x ::> x ::> x ::> Ok
@@ -588,6 +669,10 @@ foreign import ccall unsafe "gemm_mod_int64_t" c_gemmML :: Z -> Tgemm Z
588 669
589-------------------------------------------------------------------------------- 670--------------------------------------------------------------------------------
590 671
672reorderAux :: (TransArray c, Storable t, Storable a1, Storable t1, Storable a) =>
673 (CInt -> Ptr a -> CInt -> Ptr t1
674 -> Trans c (CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt))
675 -> Vector t1 -> c -> Vector t -> Vector a1
591reorderAux f s d v = unsafePerformIO $ do 676reorderAux f s d v = unsafePerformIO $ do
592 k <- createVector (dim s) 677 k <- createVector (dim s)
593 r <- createVector (dim v) 678 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 @@
13{-# LANGUAGE TypeFamilies #-} 13{-# LANGUAGE TypeFamilies #-}
14{-# LANGUAGE TypeOperators #-} 14{-# LANGUAGE TypeOperators #-}
15 15
16{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
17{-# OPTIONS_GHC -fno-warn-missing-methods #-}
18
16{- | 19{- |
17Module : Internal.Modular 20Module : Internal.Modular
18Copyright : (c) Alberto Ruiz 2015 21Copyright : (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 @@
5{-# LANGUAGE FunctionalDependencies #-} 5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-} 6{-# LANGUAGE UndecidableInstances #-}
7 7
8{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
9
8----------------------------------------------------------------------------- 10-----------------------------------------------------------------------------
9-- | 11-- |
10-- Module : Data.Packed.Internal.Numeric 12-- Module : Data.Packed.Internal.Numeric
@@ -788,13 +790,7 @@ type instance RealOf (Complex Float) = Float
788type instance RealOf I = I 790type instance RealOf I = I
789type instance RealOf Z = Z 791type instance RealOf Z = Z
790 792
791type family ComplexOf x 793type ComplexOf x = Complex (RealOf x)
792
793type instance ComplexOf Double = Complex Double
794type instance ComplexOf (Complex Double) = Complex Double
795
796type instance ComplexOf Float = Complex Float
797type instance ComplexOf (Complex Float) = Complex Float
798 794
799type family SingleOf x 795type family SingleOf x
800 796
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)
81unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x 81unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
82 82
83{-# INLINE safeIndexV #-} 83{-# INLINE safeIndexV #-}
84safeIndexV :: Storable t2
85 => (STVector s t2 -> Int -> t) -> STVector t1 t2 -> Int -> t
84safeIndexV f (STVector v) k 86safeIndexV f (STVector v) k
85 | k < 0 || k>= dim v = error $ "out of range error in vector (dim=" 87 | k < 0 || k>= dim v = error $ "out of range error in vector (dim="
86 ++show (dim v)++", pos="++show k++")" 88 ++show (dim v)++", pos="++show k++")"
@@ -150,9 +152,12 @@ unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
150freezeMatrix :: (Element t) => STMatrix s t -> ST s (Matrix t) 152freezeMatrix :: (Element t) => STMatrix s t -> ST s (Matrix t)
151freezeMatrix m = liftSTMatrix id m 153freezeMatrix m = liftSTMatrix id m
152 154
155cloneMatrix :: Element t => Matrix t -> IO (Matrix t)
153cloneMatrix m = copy (orderOf m) m 156cloneMatrix m = copy (orderOf m) m
154 157
155{-# INLINE safeIndexM #-} 158{-# INLINE safeIndexM #-}
159safeIndexM :: (STMatrix s t2 -> Int -> Int -> t)
160 -> STMatrix t1 t2 -> Int -> Int -> t
156safeIndexM f (STMatrix m) r c 161safeIndexM f (STMatrix m) r c
157 | r<0 || r>=rows m || 162 | r<0 || r>=rows m ||
158 c<0 || c>=cols m = error $ "out of range error in matrix (size=" 163 c<0 || c>=cols m = error $ "out of range error in matrix (size="
@@ -184,6 +189,7 @@ data ColRange = AllCols
184 | Col Int 189 | Col Int
185 | FromCol Int 190 | FromCol Int
186 191
192getColRange :: Int -> ColRange -> (Int, Int)
187getColRange c AllCols = (0,c-1) 193getColRange c AllCols = (0,c-1)
188getColRange c (ColRange a b) = (a `mod` c, b `mod` c) 194getColRange c (ColRange a b) = (a `mod` c, b `mod` c)
189getColRange c (Col a) = (a `mod` c, a `mod` c) 195getColRange c (Col a) = (a `mod` c, a `mod` c)
@@ -194,6 +200,7 @@ data RowRange = AllRows
194 | Row Int 200 | Row Int
195 | FromRow Int 201 | FromRow Int
196 202
203getRowRange :: Int -> RowRange -> (Int, Int)
197getRowRange r AllRows = (0,r-1) 204getRowRange r AllRows = (0,r-1)
198getRowRange r (RowRange a b) = (a `mod` r, b `mod` r) 205getRowRange r (RowRange a b) = (a `mod` r, b `mod` r)
199getRowRange r (Row a) = (a `mod` r, a `mod` r) 206getRowRange 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
223 i2' = i2 `mod` (rows m) 230 i2' = i2 `mod` (rows m)
224 231
225 232
233extractMatrix :: Element a => STMatrix t a -> RowRange -> ColRange -> ST s (Matrix a)
226extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[i1,i2]) 0 (idxs[j1,j2])) 234extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[i1,i2]) 0 (idxs[j1,j2]))
227 where 235 where
228 (i1,i2) = getRowRange (rows m) rr 236 (i1,i2) = getRowRange (rows m) rr
@@ -231,6 +239,7 @@ extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[
231-- | r0 c0 height width 239-- | r0 c0 height width
232data Slice s t = Slice (STMatrix s t) Int Int Int Int 240data Slice s t = Slice (STMatrix s t) Int Int Int Int
233 241
242slice :: Element a => Slice t a -> Matrix a
234slice (Slice (STMatrix m) r0 c0 nr nc) = subMatrix (r0,c0) (nr,nc) m 243slice (Slice (STMatrix m) r0 c0 nr nc) = subMatrix (r0,c0) (nr,nc) m
235 244
236gemmm :: Element t => t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s () 245gemmm :: 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
238 where 247 where
239 res = unsafeIOToST (gemm v a b r) 248 res = unsafeIOToST (gemm v a b r)
240 v = fromList [alpha,beta] 249 v = fromList [alpha,beta]
241 250
242 251
243mutable :: Element t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u) 252mutable :: Element t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u)
244mutable f a = runST $ do 253mutable f a = runST $ do
@@ -246,4 +255,3 @@ mutable f a = runST $ do
246 info <- f (rows a, cols a) x 255 info <- f (rows a, cols a) x
247 r <- unsafeFreezeMatrix x 256 r <- unsafeFreezeMatrix x
248 return (r,info) 257 return (r,info)
249
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 @@
2{-# LANGUAGE MultiParamTypeClasses #-} 2{-# LANGUAGE MultiParamTypeClasses #-}
3{-# LANGUAGE FlexibleInstances #-} 3{-# LANGUAGE FlexibleInstances #-}
4 4
5{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
6
5module Internal.Sparse( 7module Internal.Sparse(
6 GMatrix(..), CSR(..), mkCSR, fromCSR, 8 GMatrix(..), CSR(..), mkCSR, fromCSR,
7 mkSparse, mkDiagR, mkDense, 9 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 @@
15{-# LANGUAGE BangPatterns #-} 15{-# LANGUAGE BangPatterns #-}
16{-# LANGUAGE DeriveGeneric #-} 16{-# LANGUAGE DeriveGeneric #-}
17 17
18{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
19{-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}
20
18{- | 21{- |
19Module : Internal.Static 22Module : Internal.Static
20Copyright : (c) Alberto Ruiz 2006-14 23Copyright : (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 @@
6{-# LANGUAGE ScopedTypeVariables #-} 6{-# LANGUAGE ScopedTypeVariables #-}
7{-# LANGUAGE ViewPatterns #-} 7{-# LANGUAGE ViewPatterns #-}
8 8
9{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
10{-# OPTIONS_GHC -fno-warn-orphans #-}
9 11
10----------------------------------------------------------------------------- 12-----------------------------------------------------------------------------
11{- | 13{- |
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 @@
1{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-} 1{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-}
2{-# LANGUAGE TypeSynonymInstances #-} 2{-# LANGUAGE TypeSynonymInstances #-}
3 3
4{-# OPTIONS_GHC -fno-warn-orphans #-}
4 5
5-- | 6-- |
6-- Module : Internal.Vector 7-- Module : Internal.Vector
@@ -40,6 +41,7 @@ import qualified Data.Vector.Storable as Vector
40import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith) 41import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith)
41 42
42import Data.Binary 43import Data.Binary
44import Data.Binary.Put
43import Control.Monad(replicateM) 45import Control.Monad(replicateM)
44import qualified Data.ByteString.Internal as BS 46import qualified Data.ByteString.Internal as BS
45import Data.Vector.Storable.Internal(updPtr) 47import Data.Vector.Storable.Internal(updPtr)
@@ -92,6 +94,7 @@ createVector n = do
92 94
93-} 95-}
94 96
97safeRead :: Storable a => Vector a -> (Ptr a -> IO c) -> c
95safeRead v = inlinePerformIO . unsafeWith v 98safeRead v = inlinePerformIO . unsafeWith v
96{-# INLINE safeRead #-} 99{-# INLINE safeRead #-}
97 100
@@ -287,11 +290,13 @@ foldVectorWithIndex f x v = unsafePerformIO $
287 go (dim v -1) x 290 go (dim v -1) x
288{-# INLINE foldVectorWithIndex #-} 291{-# INLINE foldVectorWithIndex #-}
289 292
293foldLoop :: (Int -> t -> t) -> t -> Int -> t
290foldLoop f s0 d = go (d - 1) s0 294foldLoop f s0 d = go (d - 1) s0
291 where 295 where
292 go 0 s = f (0::Int) s 296 go 0 s = f (0::Int) s
293 go !j !s = go (j - 1) (f j s) 297 go !j !s = go (j - 1) (f j s)
294 298
299foldVectorG :: Storable t1 => (Int -> (Int -> t1) -> t -> t) -> t -> Vector t1 -> t
295foldVectorG f s0 v = foldLoop g s0 (dim v) 300foldVectorG f s0 v = foldLoop g s0 (dim v)
296 where g !k !s = f k (safeRead v . flip peekElemOff) s 301 where g !k !s = f k (safeRead v . flip peekElemOff) s
297 {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) 302 {-# 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
394 m = d `mod` chunk 399 m = d `mod` chunk
395 in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk) 400 in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk)
396 401
402putVector :: (Storable t, Binary t) => Vector t -> Data.Binary.Put.PutM ()
397putVector v = mapM_ put $! toList v 403putVector v = mapM_ put $! toList v
398 404
405getVector :: (Storable a, Binary a) => Int -> Get (Vector a)
399getVector d = do 406getVector d = do
400 xs <- replicateM d get 407 xs <- replicateM d get
401 return $! fromList xs 408 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)
28import Control.Monad(when) 28import Control.Monad(when)
29 29
30infixr 1 # 30infixr 1 #
31(#) :: TransArray c => c -> (b -> IO r) -> TransRaw c b -> IO r
31a # b = applyRaw a b 32a # b = applyRaw a b
32{-# INLINE (#) #-} 33{-# INLINE (#) #-}
33 34
35(#!) :: (TransArray c, TransArray c1) => c1 -> c -> TransRaw c1 (TransRaw c (IO r)) -> IO r
34a #! b = a # b # id 36a #! b = a # b # id
35{-# INLINE (#!) #-} 37{-# INLINE (#!) #-}
36 38
39fromei :: Enum a => a -> CInt
37fromei x = fromIntegral (fromEnum x) :: CInt 40fromei x = fromIntegral (fromEnum x) :: CInt
38 41
39data FunCodeV = Sin 42data FunCodeV = Sin
@@ -100,10 +103,20 @@ sumQ = sumg c_sumQ
100sumC :: Vector (Complex Double) -> Complex Double 103sumC :: Vector (Complex Double) -> Complex Double
101sumC = sumg c_sumC 104sumC = sumg c_sumC
102 105
106sumI :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr I -> I :> Ok)
107 , TransArray c
108 , Storable a
109 )
110 => I -> c -> a
103sumI m = sumg (c_sumI m) 111sumI m = sumg (c_sumI m)
104 112
113sumL :: ( TransRaw c (CInt -> Ptr a -> IO CInt) ~ (CInt -> Ptr Z -> Z :> Ok)
114 , TransArray c
115 , Storable a
116 ) => Z -> c -> a
105sumL m = sumg (c_sumL m) 117sumL m = sumg (c_sumL m)
106 118
119sumg :: (TransArray c, Storable a) => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a
107sumg f x = unsafePerformIO $ do 120sumg f x = unsafePerformIO $ do
108 r <- createVector 1 121 r <- createVector 1
109 (x #! r) f #| "sum" 122 (x #! r) f #| "sum"
@@ -140,6 +153,8 @@ prodI = prodg . c_prodI
140prodL :: Z-> Vector Z -> Z 153prodL :: Z-> Vector Z -> Z
141prodL = prodg . c_prodL 154prodL = prodg . c_prodL
142 155
156prodg :: (TransArray c, Storable a)
157 => TransRaw c (CInt -> Ptr a -> IO CInt) -> c -> a
143prodg f x = unsafePerformIO $ do 158prodg f x = unsafePerformIO $ do
144 r <- createVector 1 159 r <- createVector 1
145 (x #! r) f #| "prod" 160 (x #! r) f #| "prod"
@@ -155,16 +170,25 @@ foreign import ccall unsafe "prodL" c_prodL :: Z -> TVV Z
155 170
156------------------------------------------------------------------ 171------------------------------------------------------------------
157 172
173toScalarAux :: (Enum a, TransArray c, Storable a1)
174 => (CInt -> TransRaw c (CInt -> Ptr a1 -> IO CInt)) -> a -> c -> a1
158toScalarAux fun code v = unsafePerformIO $ do 175toScalarAux fun code v = unsafePerformIO $ do
159 r <- createVector 1 176 r <- createVector 1
160 (v #! r) (fun (fromei code)) #|"toScalarAux" 177 (v #! r) (fun (fromei code)) #|"toScalarAux"
161 return (r @> 0) 178 return (r @> 0)
162 179
180
181vectorMapAux :: (Enum a, Storable t, Storable a1)
182 => (CInt -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)
183 -> a -> Vector t -> Vector a1
163vectorMapAux fun code v = unsafePerformIO $ do 184vectorMapAux fun code v = unsafePerformIO $ do
164 r <- createVector (dim v) 185 r <- createVector (dim v)
165 (v #! r) (fun (fromei code)) #|"vectorMapAux" 186 (v #! r) (fun (fromei code)) #|"vectorMapAux"
166 return r 187 return r
167 188
189vectorMapValAux :: (Enum a, Storable a2, Storable t, Storable a1)
190 => (CInt -> Ptr a2 -> CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt)
191 -> a -> a2 -> Vector t -> Vector a1
168vectorMapValAux fun code val v = unsafePerformIO $ do 192vectorMapValAux fun code val v = unsafePerformIO $ do
169 r <- createVector (dim v) 193 r <- createVector (dim v)
170 pval <- newArray [val] 194 pval <- newArray [val]
@@ -172,6 +196,9 @@ vectorMapValAux fun code val v = unsafePerformIO $ do
172 free pval 196 free pval
173 return r 197 return r
174 198
199vectorZipAux :: (Enum a, TransArray c, Storable t, Storable a1)
200 => (CInt -> CInt -> Ptr t -> TransRaw c (CInt -> Ptr a1 -> IO CInt))
201 -> a -> Vector t -> c -> Vector a1
175vectorZipAux fun code u v = unsafePerformIO $ do 202vectorZipAux fun code u v = unsafePerformIO $ do
176 r <- createVector (dim u) 203 r <- createVector (dim u)
177 (u # v #! r) (fun (fromei code)) #|"vectorZipAux" 204 (u # v #! r) (fun (fromei code)) #|"vectorZipAux"
@@ -378,6 +405,7 @@ foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> D
378 405
379-------------------------------------------------------------------------------- 406--------------------------------------------------------------------------------
380 407
408roundVector :: Vector Double -> Vector Double
381roundVector v = unsafePerformIO $ do 409roundVector v = unsafePerformIO $ do
382 r <- createVector (dim v) 410 r <- createVector (dim v)
383 (v #! r) c_round_vector #|"roundVector" 411 (v #! r) c_round_vector #|"roundVector"
@@ -433,6 +461,8 @@ long2intV :: Vector Z -> Vector I
433long2intV = tog c_long2int 461long2intV = tog c_long2int
434 462
435 463
464tog :: (Storable t, Storable a)
465 => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
436tog f v = unsafePerformIO $ do 466tog f v = unsafePerformIO $ do
437 r <- createVector (dim v) 467 r <- createVector (dim v)
438 (v #! r) f #|"tog" 468 (v #! r) f #|"tog"
@@ -452,6 +482,8 @@ foreign import ccall unsafe "long2int" c_long2int :: Z :> I :> Ok
452 482
453--------------------------------------------------------------- 483---------------------------------------------------------------
454 484
485stepg :: (Storable t, Storable a)
486 => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
455stepg f v = unsafePerformIO $ do 487stepg f v = unsafePerformIO $ do
456 r <- createVector (dim v) 488 r <- createVector (dim v)
457 (v #! r) f #|"step" 489 (v #! r) f #|"step"
@@ -477,6 +509,8 @@ foreign import ccall unsafe "stepL" c_stepL :: TVV Z
477 509
478-------------------------------------------------------------------------------- 510--------------------------------------------------------------------------------
479 511
512conjugateAux :: (Storable t, Storable a)
513 => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
480conjugateAux fun x = unsafePerformIO $ do 514conjugateAux fun x = unsafePerformIO $ do
481 v <- createVector (dim x) 515 v <- createVector (dim x)
482 (x #! v) fun #|"conjugateAux" 516 (x #! v) fun #|"conjugateAux"
@@ -502,6 +536,8 @@ cloneVector v = do
502 536
503-------------------------------------------------------------------------------- 537--------------------------------------------------------------------------------
504 538
539constantAux :: (Storable a1, Storable a)
540 => (Ptr a1 -> CInt -> Ptr a -> IO CInt) -> a1 -> Int -> Vector a
505constantAux fun x n = unsafePerformIO $ do 541constantAux fun x n = unsafePerformIO $ do
506 v <- createVector n 542 v <- createVector n
507 px <- newArray [x] 543 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2{-# LANGUAGE FlexibleContexts #-} 2{-# LANGUAGE FlexibleContexts #-}
3 3
4{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
5
4----------------------------------------------------------------------------- 6-----------------------------------------------------------------------------
5{- | 7{- |
6Module : Numeric.LinearAlgebra 8Module : 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 <·>
28(<·>) :: Numeric t => Vector t -> Vector t -> t 28(<·>) :: Numeric t => Vector t -> Vector t -> t
29(<·>) = dot 29(<·>) = dot
30 30
31app :: Numeric t => Matrix t -> Vector t -> Vector t
31app m v = m #> v 32app m v = m #> v
32 33
34mul :: Numeric t => Matrix t -> Matrix t -> Matrix t
33mul a b = a <> b 35mul a b = a <> b
34 36
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 @@
14{-# LANGUAGE GADTs #-} 14{-# LANGUAGE GADTs #-}
15{-# LANGUAGE TypeFamilies #-} 15{-# LANGUAGE TypeFamilies #-}
16 16
17{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
18{-# OPTIONS_GHC -fno-warn-orphans #-}
17 19
18{- | 20{- |
19Module : Numeric.LinearAlgebra.Static 21Module : 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 @@
4{-# LANGUAGE UndecidableInstances #-} 4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-} 5{-# LANGUAGE MultiParamTypeClasses #-}
6 6
7{-# OPTIONS_GHC -fno-warn-orphans #-}
8
7----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
8-- | 10-- |
9-- Module : Numeric.Matrix 11-- Module : Numeric.Matrix
@@ -35,6 +37,7 @@ import Data.List(partition)
35import qualified Data.Foldable as F 37import qualified Data.Foldable as F
36import qualified Data.Semigroup as S 38import qualified Data.Semigroup as S
37import Internal.Chain 39import Internal.Chain
40import Foreign.Storable(Storable)
38 41
39 42
40------------------------------------------------------------------- 43-------------------------------------------------------------------
@@ -80,8 +83,16 @@ instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matri
80 83
81-------------------------------------------------------------------------------- 84--------------------------------------------------------------------------------
82 85
86isScalar :: Matrix t -> Bool
83isScalar m = rows m == 1 && cols m == 1 87isScalar m = rows m == 1 && cols m == 1
84 88
89adaptScalarM :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2)
90 => (t1 -> Matrix t2 -> t)
91 -> (Matrix t1 -> Matrix t2 -> t)
92 -> (Matrix t1 -> t2 -> t)
93 -> Matrix t1
94 -> Matrix t2
95 -> t
85adaptScalarM f1 f2 f3 x y 96adaptScalarM f1 f2 f3 x y
86 | isScalar x = f1 (x @@>(0,0) ) y 97 | isScalar x = f1 (x @@>(0,0) ) y
87 | isScalar y = f3 x (y @@>(0,0) ) 98 | 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
96 where 107 where
97 mempty = 1 108 mempty = 1
98 mappend = adaptScalarM scale mXm (flip scale) 109 mappend = adaptScalarM scale mXm (flip scale)
99 110
100 mconcat xs = work (partition isScalar xs) 111 mconcat xs = work (partition isScalar xs)
101 where 112 where
102 work (ss,[]) = product ss 113 work (ss,[]) = product ss
@@ -106,4 +117,3 @@ instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matr
106 | otherwise = scale x00 m 117 | otherwise = scale x00 m
107 where 118 where
108 x00 = x @@> (0,0) 119 x00 = x @@> (0,0)
109
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 @@
3{-# LANGUAGE FlexibleInstances #-} 3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE UndecidableInstances #-} 4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-} 5{-# LANGUAGE MultiParamTypeClasses #-}
6
7{-# OPTIONS_GHC -fno-warn-orphans #-}
8
6----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
7-- | 10-- |
8-- Module : Numeric.Vector 11-- Module : Numeric.Vector
@@ -14,7 +17,7 @@
14-- 17--
15-- Provides instances of standard classes 'Show', 'Read', 'Eq', 18-- Provides instances of standard classes 'Show', 'Read', 'Eq',
16-- 'Num', 'Fractional', and 'Floating' for 'Vector'. 19-- 'Num', 'Fractional', and 'Floating' for 'Vector'.
17-- 20--
18----------------------------------------------------------------------------- 21-----------------------------------------------------------------------------
19 22
20module Numeric.Vector () where 23module Numeric.Vector () where
@@ -23,9 +26,17 @@ import Internal.Vectorized
23import Internal.Vector 26import Internal.Vector
24import Internal.Numeric 27import Internal.Numeric
25import Internal.Conversion 28import Internal.Conversion
29import Foreign.Storable(Storable)
26 30
27------------------------------------------------------------------- 31-------------------------------------------------------------------
28 32
33adaptScalar :: (Foreign.Storable.Storable t1, Foreign.Storable.Storable t2)
34 => (t1 -> Vector t2 -> t)
35 -> (Vector t1 -> Vector t2 -> t)
36 -> (Vector t1 -> t2 -> t)
37 -> Vector t1
38 -> Vector t2
39 -> t
29adaptScalar f1 f2 f3 x y 40adaptScalar f1 f2 f3 x y
30 | dim x == 1 = f1 (x@>0) y 41 | dim x == 1 = f1 (x@>0) y
31 | dim y == 1 = f3 x (y@>0) 42 | dim y == 1 = f3 x (y@>0)
@@ -172,4 +183,3 @@ instance Floating (Vector (Complex Float)) where
172 sqrt = vectorMapQ Sqrt 183 sqrt = vectorMapQ Sqrt
173 (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) 184 (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
174 pi = fromList [pi] 185 pi = fromList [pi]
175
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 @@
1Name: hmatrix-glpk 1Name: hmatrix-glpk
2Version: 0.6.0.0 2Version: 0.19.0.0
3License: GPL 3License: GPL-3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
6Maintainer: Alberto Ruiz <aruiz@um.es> 6Maintainer: 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 @@
2-------- 2--------
3 3
4 * Added interpolation modules 4 * Added interpolation modules
5 5
6 * Added simulated annealing module 6 * Added simulated annealing module
7 7
8 * Added odeSolveVWith 8 * Added odeSolveVWith
@@ -11,4 +11,3 @@
11-------- 11--------
12 12
13 * The modules Numeric.GSL.* have been moved from hmatrix to the new package hmatrix-gsl. 13 * The modules Numeric.GSL.* have been moved from hmatrix to the new package hmatrix-gsl.
14
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 @@
1Name: hmatrix-gsl 1Name: hmatrix-gsl
2Version: 0.18.0.1 2Version: 0.19.0.0
3License: GPL
4License-file: LICENSE
5Author: Alberto Ruiz
6Maintainer: Alberto Ruiz <aruiz@um.es>
7Stability: provisional
8Homepage: https://github.com/albertoruiz/hmatrix
9Synopsis: Numerical computation 3Synopsis: Numerical computation
10Description: Purely functional interface to selected numerical computations, 4Description: Purely functional interface to selected numerical computations,
11 internally implemented using GSL. 5 internally implemented using GSL.
12 6Homepage: https://github.com/albertoruiz/hmatrix
7license: GPL-3
8license-file: LICENSE
9Author: Alberto Ruiz
10Maintainer: Alberto Ruiz <aruiz@um.es>
11Stability: provisional
13Category: Math 12Category: Math
14tested-with: GHC ==7.8
15
16cabal-version: >=1.8
17
18build-type: Simple 13build-type: Simple
14cabal-version: >=1.18
15
19 16
20extra-source-files: src/Numeric/GSL/gsl-ode.c 17extra-source-files: src/Numeric/GSL/gsl-ode.c
21 18
@@ -33,9 +30,6 @@ library
33 Build-Depends: base<5, hmatrix>=0.18, array, vector, 30 Build-Depends: base<5, hmatrix>=0.18, array, vector,
34 process, random 31 process, random
35 32
36
37 Extensions: ForeignFunctionInterface
38
39 hs-source-dirs: src 33 hs-source-dirs: src
40 Exposed-modules: Numeric.GSL.Differentiation, 34 Exposed-modules: Numeric.GSL.Differentiation,
41 Numeric.GSL.Integration, 35 Numeric.GSL.Integration,
@@ -98,6 +92,8 @@ library
98 else 92 else
99 pkgconfig-depends: gsl 93 pkgconfig-depends: gsl
100 94
95 default-language: Haskell2010
96
101 97
102source-repository head 98source-repository head
103 type: git 99 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 @@
1{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
2
1----------------------------------------------------------------------------- 3-----------------------------------------------------------------------------
2-- | 4-- |
3-- Module : Graphics.Plot 5-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2{-# LANGUAGE FlexibleContexts #-} 2{-# LANGUAGE FlexibleContexts #-}
3 3
4{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
5
4{- | 6{- |
5Module : Numeric.GSL.Fitting 7Module : Numeric.GSL.Fitting
6Copyright : (c) Alberto Ruiz 2010 8Copyright : (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 @@
1{-# LANGUAGE TypeFamilies #-} 1{-# LANGUAGE TypeFamilies #-}
2 2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4
3{- | 5{- |
4Module : Numeric.GSL.Fourier 6Module : Numeric.GSL.Fourier
5Copyright : (c) Alberto Ruiz 2006 7Copyright : (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 @@
1{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Integration 4Module : Numeric.GSL.Integration
3Copyright : (c) Alberto Ruiz 2006 5Copyright : (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 @@
1{-# LANGUAGE FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2 2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
5
3-- | 6-- |
4-- Module : Numeric.GSL.Internal 7-- Module : Numeric.GSL.Internal
5-- Copyright : (c) Alberto Ruiz 2009 8-- Copyright : (c) Alberto Ruiz 2009
@@ -128,8 +131,7 @@ type TVM = TV (TM Res)
128ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 131ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
129 132
130vec x f = unsafeWith x $ \p -> do 133vec x f = unsafeWith x $ \p -> do
131 let v g = do 134 let v g = g (fi $ V.length x) p
132 g (fi $ V.length x) p
133 f v 135 f v
134{-# INLINE vec #-} 136{-# INLINE vec #-}
135 137
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 @@
1{-# LANGUAGE MagicHash, UnboxedTuples #-} 1{-# LANGUAGE MagicHash, UnboxedTuples #-}
2 2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4
3{- | 5{- |
4Module : Numeric.GSL.Interpolation 6Module : Numeric.GSL.Interpolation
5Copyright : (c) Matthew Peddie 2015 7Copyright : (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 @@
1{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
2{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
3
1----------------------------------------------------------------------------- 4-----------------------------------------------------------------------------
2-- | 5-- |
3-- Module : Numeric.GSL.LinearAlgebra 6-- 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 @@
1{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
2
1{-# LANGUAGE FlexibleContexts #-} 3{-# LANGUAGE FlexibleContexts #-}
2 4
3 5
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 @@
1{-# LANGUAGE FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2 2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
3 5
4{- | 6{- |
5Module : Numeric.GSL.ODE 7Module : 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 @@
1{-# LANGUAGE FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2 2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4
3{- | 5{- |
4Module : Numeric.GSL.Root 6Module : Numeric.GSL.Root
5Copyright : (c) Alberto Ruiz 2009 7Copyright : (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 @@
1{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
2{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
3
1----------------------------------------------------------------------------- 4-----------------------------------------------------------------------------
2-- | 5-- |
3-- Module : Numeric.GSL.Vector 6-- 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 @@
1Name: hmatrix-sparse 1Name: hmatrix-sparse
2Version: 0.1.0 2Version: 0.19.0.0
3License: BSD3 3License: BSD3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
@@ -19,7 +19,7 @@ build-type: Simple
19 19
20 20
21library 21library
22 Build-Depends: base, hmatrix>=0.16 22 Build-Depends: base<5, hmatrix>=0.16
23 23
24 hs-source-dirs: src 24 hs-source-dirs: src
25 25
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 @@
1Name: hmatrix-special 1Name: hmatrix-special
2Version: 0.4.0.1 2Version: 0.19.0.0
3License: GPL 3License: GPL-3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
6Maintainer: Alberto Ruiz <aruiz@um.es> 6Maintainer: 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Bessel 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Coulomb 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Coupling 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Exp 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Gamma 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Gegenbauer 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
4
2------------------------------------------------------------ 5------------------------------------------------------------
3-- | 6-- |
4-- Module : Numeric.GSL.Special.Legendre 7-- 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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2
3{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
4{-# OPTIONS_GHC -fno-warn-unused-top-binds #-}
5
2------------------------------------------------------------ 6------------------------------------------------------------
3-- | 7-- |
4-- Module : Numeric.GSL.Special.Trig 8-- 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 @@
1# Revision history for hmatrix-sundials
2
3## 0.1.0.0 -- 2018-04-21
4
5* 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 @@
1Copyright (c) 2018, Dominic Steinitz, Novadiscovery
2
3All rights reserved.
4
5Redistribution and use in source and binary forms, with or without
6modification, are permitted provided that the following conditions are met:
7
8 * Redistributions of source code must retain the above copyright
9 notice, this list of conditions and the following disclaimer.
10
11 * Redistributions in binary form must reproduce the above
12 copyright notice, this list of conditions and the following
13 disclaimer in the documentation and/or other materials provided
14 with the distribution.
15
16 * Neither the name of Dominic Steinitz nor the names of other
17 contributors may be used to endorse or promote products derived
18 from this software without specific prior written permission.
19
20THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30OF 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 @@
1Currently only an interface to the Runge-Kutta methods:
2[ARKode](https://computation.llnl.gov/projects/sundials/arkode)
3
4The interface is almost certainly going to change. Sundials gives a
5rich set of "combinators" for controlling the solution of your problem
6and reporting on how it performed. The idea is to initially mimic
7hmatrix-gsl and add extra, richer functions but ultimately upgrade the
8whole 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 @@
1import Distribution.Simple
2main = 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 @@
1name: hmatrix-sundials
2version: 0.19.0.0
3synopsis: hmatrix interface to sundials
4description: An interface to the solving suite SUNDIALS. Currently, it
5 mimics the solving interface in hmstrix-gsl but
6 provides more diagnostic information and the
7 Butcher Tableaux (for Runge-Kutta methods).
8homepage: https://github.com/idontgetoutmuch/hmatrix/tree/sundials
9license: BSD3
10license-file: LICENSE
11author: Dominic Steinitz
12maintainer: dominic@steinitz.org
13copyright: Dominic Steinitz 2018, Novadiscovery 2018
14category: Math
15build-type: Simple
16extra-source-files: ChangeLog.md, README.md, diagrams/*.png
17extra-doc-files: diagrams/*.png
18cabal-version: >=1.18
19
20
21library
22 build-depends: base >=4.10 && <4.11,
23 inline-c >=0.6 && <0.7,
24 vector >=0.12 && <0.13,
25 template-haskell >=2.12 && <2.13,
26 containers >=0.5 && <0.6,
27 hmatrix>=0.18
28 extra-libraries: sundials_arkode
29 other-extensions: QuasiQuotes
30 hs-source-dirs: src
31 exposed-modules: Numeric.Sundials.ARKode.ODE
32 other-modules: Types,
33 Arkode
34 c-sources: src/helpers.c src/helpers.h
35 default-language: Haskell2010
36
37test-suite hmatrix-sundials-testsuite
38 type: exitcode-stdio-1.0
39 main-is: Main.hs
40 other-modules: Types,
41 Numeric.Sundials.ARKode.ODE,
42 Arkode
43 build-depends: base >=4.10 && <4.11,
44 inline-c >=0.6 && <0.7,
45 vector >=0.12 && <0.13,
46 template-haskell >=2.12 && <2.13,
47 containers >=0.5 && <0.6,
48 hmatrix>=0.18,
49 plots,
50 diagrams-lib,
51 diagrams-rasterific,
52 lens,
53 hspec
54 hs-source-dirs: src
55 extra-libraries: sundials_arkode
56 c-sources: src/helpers.c src/helpers.h
57 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 @@
1module Arkode where
2
3import Foreign
4import Foreign.C.Types
5
6
7#include <stdio.h>
8#include <sundials/sundials_nvector.h>
9#include <sundials/sundials_matrix.h>
10#include <nvector/nvector_serial.h>
11#include <sunmatrix/sunmatrix_dense.h>
12#include <arkode/arkode.h>
13
14
15#def typedef struct _generic_N_Vector SunVector;
16#def typedef struct _N_VectorContent_Serial SunContent;
17
18#def typedef struct _generic_SUNMatrix SunMatrix;
19#def typedef struct _SUNMatrixContent_Dense SunMatrixContent;
20
21getContentMatrixPtr :: Storable a => Ptr b -> IO a
22getContentMatrixPtr ptr = (#peek SunMatrix, content) ptr
23
24getNRows :: Ptr b -> IO CInt
25getNRows ptr = (#peek SunMatrixContent, M) ptr
26putNRows :: CInt -> Ptr b -> IO ()
27putNRows nr ptr = (#poke SunMatrixContent, M) ptr nr
28
29getNCols :: Ptr b -> IO CInt
30getNCols ptr = (#peek SunMatrixContent, N) ptr
31putNCols :: CInt -> Ptr b -> IO ()
32putNCols nc ptr = (#poke SunMatrixContent, N) ptr nc
33
34getMatrixData :: Storable a => Ptr b -> IO a
35getMatrixData ptr = (#peek SunMatrixContent, data) ptr
36
37getContentPtr :: Storable a => Ptr b -> IO a
38getContentPtr ptr = (#peek SunVector, content) ptr
39
40getData :: Storable a => Ptr b -> IO a
41getData ptr = (#peek SunContent, data) ptr
42
43arkSMax :: Int
44arkSMax = #const ARK_S_MAX
45
46mIN_DIRK_NUM, mAX_DIRK_NUM :: Int
47mIN_DIRK_NUM = #const MIN_DIRK_NUM
48mAX_DIRK_NUM = #const MAX_DIRK_NUM
49
50-- FIXME: We could just use inline-c instead
51
52-- Butcher table accessors -- implicit
53sDIRK_2_1_2 :: Int
54sDIRK_2_1_2 = #const SDIRK_2_1_2
55bILLINGTON_3_3_2 :: Int
56bILLINGTON_3_3_2 = #const BILLINGTON_3_3_2
57tRBDF2_3_3_2 :: Int
58tRBDF2_3_3_2 = #const TRBDF2_3_3_2
59kVAERNO_4_2_3 :: Int
60kVAERNO_4_2_3 = #const KVAERNO_4_2_3
61aRK324L2SA_DIRK_4_2_3 :: Int
62aRK324L2SA_DIRK_4_2_3 = #const ARK324L2SA_DIRK_4_2_3
63cASH_5_2_4 :: Int
64cASH_5_2_4 = #const CASH_5_2_4
65cASH_5_3_4 :: Int
66cASH_5_3_4 = #const CASH_5_3_4
67sDIRK_5_3_4 :: Int
68sDIRK_5_3_4 = #const SDIRK_5_3_4
69kVAERNO_5_3_4 :: Int
70kVAERNO_5_3_4 = #const KVAERNO_5_3_4
71aRK436L2SA_DIRK_6_3_4 :: Int
72aRK436L2SA_DIRK_6_3_4 = #const ARK436L2SA_DIRK_6_3_4
73kVAERNO_7_4_5 :: Int
74kVAERNO_7_4_5 = #const KVAERNO_7_4_5
75aRK548L2SA_DIRK_8_4_5 :: Int
76aRK548L2SA_DIRK_8_4_5 = #const ARK548L2SA_DIRK_8_4_5
77
78-- #define DEFAULT_DIRK_2 SDIRK_2_1_2
79-- #define DEFAULT_DIRK_3 ARK324L2SA_DIRK_4_2_3
80-- #define DEFAULT_DIRK_4 SDIRK_5_3_4
81-- #define DEFAULT_DIRK_5 ARK548L2SA_DIRK_8_4_5
82
83-- Butcher table accessors -- explicit
84hEUN_EULER_2_1_2 :: Int
85hEUN_EULER_2_1_2 = #const HEUN_EULER_2_1_2
86bOGACKI_SHAMPINE_4_2_3 :: Int
87bOGACKI_SHAMPINE_4_2_3 = #const BOGACKI_SHAMPINE_4_2_3
88aRK324L2SA_ERK_4_2_3 :: Int
89aRK324L2SA_ERK_4_2_3 = #const ARK324L2SA_ERK_4_2_3
90zONNEVELD_5_3_4 :: Int
91zONNEVELD_5_3_4 = #const ZONNEVELD_5_3_4
92aRK436L2SA_ERK_6_3_4 :: Int
93aRK436L2SA_ERK_6_3_4 = #const ARK436L2SA_ERK_6_3_4
94sAYFY_ABURUB_6_3_4 :: Int
95sAYFY_ABURUB_6_3_4 = #const SAYFY_ABURUB_6_3_4
96cASH_KARP_6_4_5 :: Int
97cASH_KARP_6_4_5 = #const CASH_KARP_6_4_5
98fEHLBERG_6_4_5 :: Int
99fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5
100dORMAND_PRINCE_7_4_5 :: Int
101dORMAND_PRINCE_7_4_5 = #const DORMAND_PRINCE_7_4_5
102aRK548L2SA_ERK_8_4_5 :: Int
103aRK548L2SA_ERK_8_4_5 = #const ARK548L2SA_ERK_8_4_5
104vERNER_8_5_6 :: Int
105vERNER_8_5_6 = #const VERNER_8_5_6
106fEHLBERG_13_7_8 :: Int
107fEHLBERG_13_7_8 = #const FEHLBERG_13_7_8
108
109-- #define DEFAULT_ERK_2 HEUN_EULER_2_1_2
110-- #define DEFAULT_ERK_3 BOGACKI_SHAMPINE_4_2_3
111-- #define DEFAULT_ERK_4 ZONNEVELD_5_3_4
112-- #define DEFAULT_ERK_5 CASH_KARP_6_4_5
113-- #define DEFAULT_ERK_6 VERNER_8_5_6
114-- #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 @@
1{-# OPTIONS_GHC -Wall #-}
2
3import Numeric.Sundials.ARKode.ODE
4import Numeric.LinearAlgebra
5
6import Plots as P
7import qualified Diagrams.Prelude as D
8import Diagrams.Backend.Rasterific
9
10import Control.Lens
11
12import Test.Hspec
13
14
15lorenz :: Double -> [Double] -> [Double]
16lorenz _t u = [ sigma * (y - x)
17 , x * (rho - z) - y
18 , x * y - beta * z
19 ]
20 where
21 rho = 28.0
22 sigma = 10.0
23 beta = 8.0 / 3.0
24 x = u !! 0
25 y = u !! 1
26 z = u !! 2
27
28_lorenzJac :: Double -> Vector Double -> Matrix Double
29_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y
30 , sigma , -1.0 , x
31 , 0.0 , (-x) , (-beta)
32 ]
33 where
34 rho = 28.0
35 sigma = 10.0
36 beta = 8.0 / 3.0
37 x = u ! 0
38 y = u ! 1
39 z = u ! 2
40
41brusselator :: Double -> [Double] -> [Double]
42brusselator _t x = [ a - (w + 1) * u + v * u * u
43 , w * u - v * u * u
44 , (b - w) / eps - w * u
45 ]
46 where
47 a = 1.0
48 b = 3.5
49 eps = 5.0e-6
50 u = x !! 0
51 v = x !! 1
52 w = x !! 2
53
54_brussJac :: Double -> Vector Double -> Matrix Double
55_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
56 , u * u , (-(u * u)) , 0.0
57 , (-u) , u , (-1.0) / eps - u
58 ]
59 where
60 y = toList x
61 u = y !! 0
62 v = y !! 1
63 w = y !! 2
64 eps = 5.0e-6
65
66stiffish :: Double -> [Double] -> [Double]
67stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
68 where
69 lamda = -100.0
70 u = v !! 0
71
72stiffishV :: Double -> Vector Double -> Vector Double
73stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
74 where
75 lamda = -100.0
76 u = v ! 0
77
78_stiffJac :: Double -> Vector Double -> Matrix Double
79_stiffJac _t _v = (1><1) [ lamda ]
80 where
81 lamda = -100.0
82
83lSaxis :: [[Double]] -> P.Axis B D.V2 Double
84lSaxis xs = P.r2Axis &~ do
85 let ts = xs!!0
86 us = xs!!1
87 vs = xs!!2
88 ws = xs!!3
89 P.linePlot' $ zip ts us
90 P.linePlot' $ zip ts vs
91 P.linePlot' $ zip ts ws
92
93kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
94kSaxis xs = P.r2Axis &~ do
95 P.linePlot' xs
96
97main :: IO ()
98main = do
99
100 let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
101 renderRasterific "diagrams/brusselator.png"
102 (D.dims2D 500.0 500.0)
103 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
104
105 let res1a = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
106 renderRasterific "diagrams/brusselatorA.png"
107 (D.dims2D 500.0 500.0)
108 (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))
109
110 let res2 = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
111 renderRasterific "diagrams/stiffish.png"
112 (D.dims2D 500.0 500.0)
113 (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2))
114
115 let res2a = odeSolveV (SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
116
117 let res2b = odeSolveV (TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])
118
119 let maxDiff = maximum $ map abs $
120 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)
121
122 hspec $ describe "Compare results" $ do
123 it "for two different RK methods" $
124 maxDiff < 1.0e-6
125
126 let res3 = odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])
127
128 renderRasterific "diagrams/lorenz.png"
129 (D.dims2D 500.0 500.0)
130 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))
131
132 renderRasterific "diagrams/lorenz1.png"
133 (D.dims2D 500.0 500.0)
134 (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2))
135
136 renderRasterific "diagrams/lorenz2.png"
137 (D.dims2D 500.0 500.0)
138 (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 @@
1{-# OPTIONS_GHC -Wall #-}
2
3{-# LANGUAGE QuasiQuotes #-}
4{-# LANGUAGE TemplateHaskell #-}
5{-# LANGUAGE MultiWayIf #-}
6{-# LANGUAGE OverloadedStrings #-}
7{-# LANGUAGE ScopedTypeVariables #-}
8{-# LANGUAGE DeriveGeneric #-}
9{-# LANGUAGE TypeOperators #-}
10{-# LANGUAGE KindSignatures #-}
11{-# LANGUAGE TypeSynonymInstances #-}
12{-# LANGUAGE FlexibleInstances #-}
13{-# LANGUAGE FlexibleContexts #-}
14
15-----------------------------------------------------------------------------
16-- |
17-- Module : Numeric.Sundials.ARKode.ODE
18-- Copyright : Dominic Steinitz 2018,
19-- Novadiscovery 2018
20-- License : BSD
21-- Maintainer : Dominic Steinitz
22-- Stability : provisional
23--
24-- Solution of ordinary differential equation (ODE) initial value problems.
25--
26-- <https://computation.llnl.gov/projects/sundials/sundials-software>
27--
28-- A simple example:
29--
30-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
31--
32-- @
33-- import Numeric.Sundials.ARKode.ODE
34-- import Numeric.LinearAlgebra
35--
36-- import Plots as P
37-- import qualified Diagrams.Prelude as D
38-- import Diagrams.Backend.Rasterific
39--
40-- brusselator :: Double -> [Double] -> [Double]
41-- brusselator _t x = [ a - (w + 1) * u + v * u * u
42-- , w * u - v * u * u
43-- , (b - w) / eps - w * u
44-- ]
45-- where
46-- a = 1.0
47-- b = 3.5
48-- eps = 5.0e-6
49-- u = x !! 0
50-- v = x !! 1
51-- w = x !! 2
52--
53-- lSaxis :: [[Double]] -> P.Axis B D.V2 Double
54-- lSaxis xs = P.r2Axis &~ do
55-- let ts = xs!!0
56-- us = xs!!1
57-- vs = xs!!2
58-- ws = xs!!3
59-- P.linePlot' $ zip ts us
60-- P.linePlot' $ zip ts vs
61-- P.linePlot' $ zip ts ws
62--
63-- main = do
64-- let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
65-- renderRasterific "diagrams/brusselator.png"
66-- (D.dims2D 500.0 500.0)
67-- (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
68-- @
69--
70-- KVAERNO_4_2_3
71--
72-- \[
73-- \begin{array}{c|cccc}
74-- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
75-- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\
76-- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
77-- 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
78-- \hline
79-- & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
80-- & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
81-- \end{array}
82-- \]
83--
84-- SDIRK_2_1_2
85--
86-- \[
87-- \begin{array}{c|cc}
88-- 1.0 & 1.0 & 0.0 \\
89-- 0.0 & -1.0 & 1.0 \\
90-- \hline
91-- & 0.5 & 0.5 \\
92-- & 1.0 & 0.0 \\
93-- \end{array}
94-- \]
95--
96-- SDIRK_5_3_4
97--
98-- \[
99-- \begin{array}{c|ccccc}
100-- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\
101-- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\
102-- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\
103-- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\
104-- 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
105-- \hline
106-- & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
107-- & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\
108-- \end{array}
109-- \]
110-----------------------------------------------------------------------------
111module Numeric.Sundials.ARKode.ODE ( odeSolve
112 , odeSolveV
113 , odeSolveVWith
114 , odeSolveVWith'
115 , ButcherTable(..)
116 , butcherTable
117 , ODEMethod(..)
118 , StepControl(..)
119 , Jacobian
120 , SundialsDiagnostics(..)
121 ) where
122
123import qualified Language.C.Inline as C
124import qualified Language.C.Inline.Unsafe as CU
125
126import Data.Monoid ((<>))
127import Data.Maybe (isJust)
128
129import Foreign.C.Types
130import Foreign.Ptr (Ptr)
131import Foreign.ForeignPtr (newForeignPtr_)
132import Foreign.Storable (Storable)
133
134import qualified Data.Vector.Storable as V
135import qualified Data.Vector.Storable.Mutable as VM
136
137import Data.Coerce (coerce)
138import System.IO.Unsafe (unsafePerformIO)
139import GHC.Generics
140
141import Numeric.LinearAlgebra.Devel (createVector)
142
143import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
144 subMatrix, rows, cols, toLists,
145 size, subVector)
146
147import qualified Types as T
148import Arkode
149import qualified Arkode as B
150
151
152C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
153
154C.include "<stdlib.h>"
155C.include "<stdio.h>"
156C.include "<math.h>"
157C.include "<arkode/arkode.h>" -- prototypes for ARKODE fcts., consts.
158C.include "<nvector/nvector_serial.h>" -- serial N_Vector types, fcts., macros
159C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
160C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
161C.include "<arkode/arkode_direct.h>" -- access to ARKDls interface
162C.include "<sundials/sundials_types.h>" -- definition of type realtype
163C.include "<sundials/sundials_math.h>"
164C.include "../../../helpers.h"
165C.include "Arkode_hsc.h"
166
167
168getDataFromContents :: Int -> Ptr T.SunVector -> IO (V.Vector CDouble)
169getDataFromContents len ptr = do
170 qtr <- B.getContentPtr ptr
171 rtr <- B.getData qtr
172 vectorFromC len rtr
173
174-- FIXME: Potentially an instance of Storable
175_getMatrixDataFromContents :: Ptr T.SunMatrix -> IO T.SunMatrix
176_getMatrixDataFromContents ptr = do
177 qtr <- B.getContentMatrixPtr ptr
178 rs <- B.getNRows qtr
179 cs <- B.getNCols qtr
180 rtr <- B.getMatrixData qtr
181 vs <- vectorFromC (fromIntegral $ rs * cs) rtr
182 return $ T.SunMatrix { T.rows = rs, T.cols = cs, T.vals = vs }
183
184putMatrixDataFromContents :: T.SunMatrix -> Ptr T.SunMatrix -> IO ()
185putMatrixDataFromContents mat ptr = do
186 let rs = T.rows mat
187 cs = T.cols mat
188 vs = T.vals mat
189 qtr <- B.getContentMatrixPtr ptr
190 B.putNRows rs qtr
191 B.putNCols cs qtr
192 rtr <- B.getMatrixData qtr
193 vectorToC vs (fromIntegral $ rs * cs) rtr
194-- FIXME: END
195
196putDataInContents :: Storable a => V.Vector a -> Int -> Ptr b -> IO ()
197putDataInContents vec len ptr = do
198 qtr <- B.getContentPtr ptr
199 rtr <- B.getData qtr
200 vectorToC vec len rtr
201
202-- Utils
203
204vectorFromC :: Storable a => Int -> Ptr a -> IO (V.Vector a)
205vectorFromC len ptr = do
206 ptr' <- newForeignPtr_ ptr
207 V.freeze $ VM.unsafeFromForeignPtr0 ptr' len
208
209vectorToC :: Storable a => V.Vector a -> Int -> Ptr a -> IO ()
210vectorToC vec len ptr = do
211 ptr' <- newForeignPtr_ ptr
212 V.copy (VM.unsafeFromForeignPtr0 ptr' len) vec
213
214data SundialsDiagnostics = SundialsDiagnostics {
215 aRKodeGetNumSteps :: Int
216 , aRKodeGetNumStepAttempts :: Int
217 , aRKodeGetNumRhsEvals_fe :: Int
218 , aRKodeGetNumRhsEvals_fi :: Int
219 , aRKodeGetNumLinSolvSetups :: Int
220 , aRKodeGetNumErrTestFails :: Int
221 , aRKodeGetNumNonlinSolvIters :: Int
222 , aRKodeGetNumNonlinSolvConvFails :: Int
223 , aRKDlsGetNumJacEvals :: Int
224 , aRKDlsGetNumRhsEvals :: Int
225 } deriving Show
226
227type Jacobian = Double -> Vector Double -> Matrix Double
228
229-- | Stepping functions
230data ODEMethod = SDIRK_2_1_2 Jacobian
231 | SDIRK_2_1_2'
232 | BILLINGTON_3_3_2 Jacobian
233 | BILLINGTON_3_3_2'
234 | TRBDF2_3_3_2 Jacobian
235 | TRBDF2_3_3_2'
236 | KVAERNO_4_2_3 Jacobian
237 | KVAERNO_4_2_3'
238 | ARK324L2SA_DIRK_4_2_3 Jacobian
239 | ARK324L2SA_DIRK_4_2_3'
240 | CASH_5_2_4 Jacobian
241 | CASH_5_2_4'
242 | CASH_5_3_4 Jacobian
243 | CASH_5_3_4'
244 | SDIRK_5_3_4 Jacobian
245 | SDIRK_5_3_4'
246 | KVAERNO_5_3_4 Jacobian
247 | KVAERNO_5_3_4'
248 | ARK436L2SA_DIRK_6_3_4 Jacobian
249 | ARK436L2SA_DIRK_6_3_4'
250 | KVAERNO_7_4_5 Jacobian
251 | KVAERNO_7_4_5'
252 | ARK548L2SA_DIRK_8_4_5 Jacobian
253 | ARK548L2SA_DIRK_8_4_5'
254 | HEUN_EULER_2_1_2 Jacobian
255 | HEUN_EULER_2_1_2'
256 | BOGACKI_SHAMPINE_4_2_3 Jacobian
257 | BOGACKI_SHAMPINE_4_2_3'
258 | ARK324L2SA_ERK_4_2_3 Jacobian
259 | ARK324L2SA_ERK_4_2_3'
260 | ZONNEVELD_5_3_4 Jacobian
261 | ZONNEVELD_5_3_4'
262 | ARK436L2SA_ERK_6_3_4 Jacobian
263 | ARK436L2SA_ERK_6_3_4'
264 | SAYFY_ABURUB_6_3_4 Jacobian
265 | SAYFY_ABURUB_6_3_4'
266 | CASH_KARP_6_4_5 Jacobian
267 | CASH_KARP_6_4_5'
268 | FEHLBERG_6_4_5 Jacobian
269 | FEHLBERG_6_4_5'
270 | DORMAND_PRINCE_7_4_5 Jacobian
271 | DORMAND_PRINCE_7_4_5'
272 | ARK548L2SA_ERK_8_4_5 Jacobian
273 | ARK548L2SA_ERK_8_4_5'
274 | VERNER_8_5_6 Jacobian
275 | VERNER_8_5_6'
276 | FEHLBERG_13_7_8 Jacobian
277 | FEHLBERG_13_7_8'
278 deriving Generic
279
280constrName :: (HasConstructor (Rep a), Generic a)=> a -> String
281constrName = genericConstrName . from
282
283class HasConstructor (f :: * -> *) where
284 genericConstrName :: f x -> String
285
286instance HasConstructor f => HasConstructor (D1 c f) where
287 genericConstrName (M1 x) = genericConstrName x
288
289instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where
290 genericConstrName (L1 l) = genericConstrName l
291 genericConstrName (R1 r) = genericConstrName r
292
293instance Constructor c => HasConstructor (C1 c f) where
294 genericConstrName x = conName x
295
296instance Show ODEMethod where
297 show x = constrName x
298
299-- FIXME: We can probably do better here with generics
300getMethod :: ODEMethod -> Int
301getMethod (SDIRK_2_1_2 _) = sDIRK_2_1_2
302getMethod (SDIRK_2_1_2') = sDIRK_2_1_2
303getMethod (BILLINGTON_3_3_2 _) = bILLINGTON_3_3_2
304getMethod (BILLINGTON_3_3_2') = bILLINGTON_3_3_2
305getMethod (TRBDF2_3_3_2 _) = tRBDF2_3_3_2
306getMethod (TRBDF2_3_3_2') = tRBDF2_3_3_2
307getMethod (KVAERNO_4_2_3 _) = kVAERNO_4_2_3
308getMethod (KVAERNO_4_2_3') = kVAERNO_4_2_3
309getMethod (ARK324L2SA_DIRK_4_2_3 _) = aRK324L2SA_DIRK_4_2_3
310getMethod (ARK324L2SA_DIRK_4_2_3') = aRK324L2SA_DIRK_4_2_3
311getMethod (CASH_5_2_4 _) = cASH_5_2_4
312getMethod (CASH_5_2_4') = cASH_5_2_4
313getMethod (CASH_5_3_4 _) = cASH_5_3_4
314getMethod (CASH_5_3_4') = cASH_5_3_4
315getMethod (SDIRK_5_3_4 _) = sDIRK_5_3_4
316getMethod (SDIRK_5_3_4') = sDIRK_5_3_4
317getMethod (KVAERNO_5_3_4 _) = kVAERNO_5_3_4
318getMethod (KVAERNO_5_3_4') = kVAERNO_5_3_4
319getMethod (ARK436L2SA_DIRK_6_3_4 _) = aRK436L2SA_DIRK_6_3_4
320getMethod (ARK436L2SA_DIRK_6_3_4') = aRK436L2SA_DIRK_6_3_4
321getMethod (KVAERNO_7_4_5 _) = kVAERNO_7_4_5
322getMethod (KVAERNO_7_4_5') = kVAERNO_7_4_5
323getMethod (ARK548L2SA_DIRK_8_4_5 _) = aRK548L2SA_DIRK_8_4_5
324getMethod (ARK548L2SA_DIRK_8_4_5') = aRK548L2SA_DIRK_8_4_5
325getMethod (HEUN_EULER_2_1_2 _) = hEUN_EULER_2_1_2
326getMethod (HEUN_EULER_2_1_2') = hEUN_EULER_2_1_2
327getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3
328getMethod (BOGACKI_SHAMPINE_4_2_3') = bOGACKI_SHAMPINE_4_2_3
329getMethod (ARK324L2SA_ERK_4_2_3 _) = aRK324L2SA_ERK_4_2_3
330getMethod (ARK324L2SA_ERK_4_2_3') = aRK324L2SA_ERK_4_2_3
331getMethod (ZONNEVELD_5_3_4 _) = zONNEVELD_5_3_4
332getMethod (ZONNEVELD_5_3_4') = zONNEVELD_5_3_4
333getMethod (ARK436L2SA_ERK_6_3_4 _) = aRK436L2SA_ERK_6_3_4
334getMethod (ARK436L2SA_ERK_6_3_4') = aRK436L2SA_ERK_6_3_4
335getMethod (SAYFY_ABURUB_6_3_4 _) = sAYFY_ABURUB_6_3_4
336getMethod (SAYFY_ABURUB_6_3_4') = sAYFY_ABURUB_6_3_4
337getMethod (CASH_KARP_6_4_5 _) = cASH_KARP_6_4_5
338getMethod (CASH_KARP_6_4_5') = cASH_KARP_6_4_5
339getMethod (FEHLBERG_6_4_5 _) = fEHLBERG_6_4_5
340getMethod (FEHLBERG_6_4_5' ) = fEHLBERG_6_4_5
341getMethod (DORMAND_PRINCE_7_4_5 _) = dORMAND_PRINCE_7_4_5
342getMethod (DORMAND_PRINCE_7_4_5') = dORMAND_PRINCE_7_4_5
343getMethod (ARK548L2SA_ERK_8_4_5 _) = aRK548L2SA_ERK_8_4_5
344getMethod (ARK548L2SA_ERK_8_4_5') = aRK548L2SA_ERK_8_4_5
345getMethod (VERNER_8_5_6 _) = vERNER_8_5_6
346getMethod (VERNER_8_5_6') = vERNER_8_5_6
347getMethod (FEHLBERG_13_7_8 _) = fEHLBERG_13_7_8
348getMethod (FEHLBERG_13_7_8') = fEHLBERG_13_7_8
349
350getJacobian :: ODEMethod -> Maybe Jacobian
351getJacobian (SDIRK_2_1_2 j) = Just j
352getJacobian (BILLINGTON_3_3_2 j) = Just j
353getJacobian (TRBDF2_3_3_2 j) = Just j
354getJacobian (KVAERNO_4_2_3 j) = Just j
355getJacobian (ARK324L2SA_DIRK_4_2_3 j) = Just j
356getJacobian (CASH_5_2_4 j) = Just j
357getJacobian (CASH_5_3_4 j) = Just j
358getJacobian (SDIRK_5_3_4 j) = Just j
359getJacobian (KVAERNO_5_3_4 j) = Just j
360getJacobian (ARK436L2SA_DIRK_6_3_4 j) = Just j
361getJacobian (KVAERNO_7_4_5 j) = Just j
362getJacobian (ARK548L2SA_DIRK_8_4_5 j) = Just j
363getJacobian (HEUN_EULER_2_1_2 j) = Just j
364getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j
365getJacobian (ARK324L2SA_ERK_4_2_3 j) = Just j
366getJacobian (ZONNEVELD_5_3_4 j) = Just j
367getJacobian (ARK436L2SA_ERK_6_3_4 j) = Just j
368getJacobian (SAYFY_ABURUB_6_3_4 j) = Just j
369getJacobian (CASH_KARP_6_4_5 j) = Just j
370getJacobian (FEHLBERG_6_4_5 j) = Just j
371getJacobian (DORMAND_PRINCE_7_4_5 j) = Just j
372getJacobian (ARK548L2SA_ERK_8_4_5 j) = Just j
373getJacobian (VERNER_8_5_6 j) = Just j
374getJacobian (FEHLBERG_13_7_8 j) = Just j
375getJacobian _ = Nothing
376
377-- | A version of 'odeSolveVWith' with reasonable default step control.
378odeSolveV
379 :: ODEMethod
380 -> Maybe Double -- ^ initial step size - by default, ARKode
381 -- estimates the initial step size to be the
382 -- solution \(h\) of the equation
383 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
384 -- \(\ddot{y}\) is an estimated value of the
385 -- second derivative of the solution at \(t_0\)
386 -> Double -- ^ absolute tolerance for the state vector
387 -> Double -- ^ relative tolerance for the state vector
388 -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
389 -> Vector Double -- ^ initial conditions
390 -> Vector Double -- ^ desired solution times
391 -> Matrix Double -- ^ solution
392odeSolveV meth hi epsAbs epsRel f y0 ts =
393 case odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts of
394 Left c -> error $ show c -- FIXME
395 -- FIXME: Can we do better than using lists?
396 Right (v, _d) -> (nR >< nC) (V.toList v)
397 where
398 us = toList ts
399 nR = length us
400 nC = size y0
401 g t x0 = coerce $ f t x0
402
403-- | A version of 'odeSolveV' with reasonable default parameters and
404-- system of equations defined using lists. FIXME: we should say
405-- something about the fact we could use the Jacobian but don't for
406-- compatibility with hmatrix-gsl.
407odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
408 -> [Double] -- ^ initial conditions
409 -> Vector Double -- ^ desired solution times
410 -> Matrix Double -- ^ solution
411odeSolve f y0 ts =
412 -- FIXME: These tolerances are different from the ones in GSL
413 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
414 Left c -> error $ show c -- FIXME
415 Right (v, _d) -> (nR >< nC) (V.toList v)
416 where
417 us = toList ts
418 nR = length us
419 nC = length y0
420 g t x0 = V.fromList $ f t (V.toList x0)
421
422odeSolveVWith' ::
423 ODEMethod
424 -> StepControl
425 -> Maybe Double -- ^ initial step size - by default, ARKode
426 -- estimates the initial step size to be the
427 -- solution \(h\) of the equation
428 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
429 -- \(\ddot{y}\) is an estimated value of the second
430 -- derivative of the solution at \(t_0\)
431 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
432 -> V.Vector Double -- ^ Initial conditions
433 -> V.Vector Double -- ^ Desired solution times
434 -> Matrix Double -- ^ Error code or solution
435odeSolveVWith' method control initStepSize f y0 tt =
436 case odeSolveVWith method control initStepSize f y0 tt of
437 Left c -> error $ show c -- FIXME
438 Right (v, _d) -> (nR >< nC) (V.toList v)
439 where
440 nR = V.length tt
441 nC = V.length y0
442
443odeSolveVWith ::
444 ODEMethod
445 -> StepControl
446 -> Maybe Double -- ^ initial step size - by default, ARKode
447 -- estimates the initial step size to be the
448 -- solution \(h\) of the equation
449 -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
450 -- \(\ddot{y}\) is an estimated value of the second
451 -- derivative of the solution at \(t_0\)
452 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
453 -> V.Vector Double -- ^ Initial conditions
454 -> V.Vector Double -- ^ Desired solution times
455 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution
456odeSolveVWith method control initStepSize f y0 tt =
457 case solveOdeC (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
458 (coerce f) (coerce y0) (coerce tt) of
459 Left c -> Left $ fromIntegral c
460 Right (v, d) -> Right (coerce v, d)
461 where
462 l = size y0
463 scise (X absTol relTol) = coerce (V.replicate l absTol, relTol)
464 scise (X' absTol relTol) = coerce (V.replicate l absTol, relTol)
465 scise (XX' absTol relTol yScale _yDotScale) = coerce (V.replicate l absTol, yScale * relTol)
466 -- FIXME; Should we check that the length of ss is correct?
467 scise (ScXX' absTol relTol yScale _yDotScale ss) = coerce (V.map (* absTol) ss, yScale * relTol)
468 jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
469 getJacobian method
470 matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
471 where
472 nr = fromIntegral $ rows m
473 nc = fromIntegral $ cols m
474 -- FIXME: efficiency
475 vs = V.fromList $ map coerce $ concat $ toLists m
476
477solveOdeC ::
478 CInt ->
479 Maybe CDouble ->
480 (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
481 (V.Vector CDouble, CDouble) ->
482 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
483 -> V.Vector CDouble -- ^ Initial conditions
484 -> V.Vector CDouble -- ^ Desired solution times
485 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
486solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO $ do
487
488 let isInitStepSize :: CInt
489 isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
490 ss :: CDouble
491 ss = case initStepSize of
492 -- It would be better to put an error message here but
493 -- inline-c seems to evaluate this even if it is never
494 -- used :(
495 Nothing -> 0.0
496 Just x -> x
497 let dim = V.length f0
498 nEq :: CLong
499 nEq = fromIntegral dim
500 nTs :: CInt
501 nTs = fromIntegral $ V.length ts
502 -- FIXME: fMut is not actually mutatated
503 fMut <- V.thaw f0
504 tMut <- V.thaw ts
505 -- FIXME: I believe this gets taken from the ghc heap and so should
506 -- be subject to garbage collection.
507 quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
508 qMatMut <- V.thaw quasiMatrixRes
509 diagnostics :: V.Vector CLong <- createVector 10 -- FIXME
510 diagMut <- V.thaw diagnostics
511 -- We need the types that sundials expects. These are tied together
512 -- in 'Types'. FIXME: The Haskell type is currently empty!
513 let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
514 funIO x y f _ptr = do
515 -- Convert the pointer we get from C (y) to a vector, and then
516 -- apply the user-supplied function.
517 fImm <- fun x <$> getDataFromContents dim y
518 -- Fill in the provided pointer with the resulting vector.
519 putDataInContents fImm dim f
520 -- FIXME: I don't understand what this comment means
521 -- Unsafe since the function will be called many times.
522 [CU.exp| int{ 0 } |]
523 let isJac :: CInt
524 isJac = fromIntegral $ fromEnum $ isJust jacH
525 jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
526 Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
527 IO CInt
528 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
529 case jacH of
530 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
531 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
532 putMatrixDataFromContents j jacS
533 -- FIXME: I don't understand what this comment means
534 -- Unsafe since the function will be called many times.
535 [CU.exp| int{ 0 } |]
536
537 res <- [C.block| int {
538 /* general problem variables */
539
540 int flag; /* reusable error-checking flag */
541 int i, j; /* reusable loop indices */
542 N_Vector y = NULL; /* empty vector for storing solution */
543 N_Vector tv = NULL; /* empty vector for storing absolute tolerances */
544 SUNMatrix A = NULL; /* empty matrix for linear solver */
545 SUNLinearSolver LS = NULL; /* empty linear solver object */
546 void *arkode_mem = NULL; /* empty ARKode memory structure */
547 realtype t;
548 long int nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;
549
550 /* general problem parameters */
551
552 realtype T0 = RCONST(($vec-ptr:(double *tMut))[0]); /* initial time */
553 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */
554
555 /* Initialize data structures */
556
557 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
558 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
559 /* Specify initial condition */
560 for (i = 0; i < NEQ; i++) {
561 NV_Ith_S(y,i) = ($vec-ptr:(double *fMut))[i];
562 };
563
564 tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
565 if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
566 /* Specify tolerances */
567 for (i = 0; i < NEQ; i++) {
568 NV_Ith_S(tv,i) = ($vec-ptr:(double *absTols))[i];
569 };
570
571 arkode_mem = ARKodeCreate(); /* Create the solver memory */
572 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
573
574 /* Call ARKodeInit to initialize the integrator memory and specify the */
575 /* right-hand side function in y'=f(t,y), the inital time T0, and */
576 /* the initial dependent variable vector y. Note: we treat the */
577 /* problem as fully implicit and set f_E to NULL and f_I to f. */
578
579 /* Here we use the C types defined in helpers.h which tie up with */
580 /* the Haskell types defined in Types */
581 if ($(int method) < MIN_DIRK_NUM) {
582 flag = ARKodeInit(arkode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y);
583 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
584 } else {
585 flag = ARKodeInit(arkode_mem, NULL, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
586 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
587 }
588
589 /* FIXME: A hack for initial testing */
590 flag = ARKodeSetMinStep(arkode_mem, 1.0e-12);
591 if (check_flag(&flag, "ARKodeSetMinStep", 1)) return 1;
592 flag = ARKodeSetMaxNumSteps(arkode_mem, 10000);
593 if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) return 1;
594
595 /* Set routines */
596 flag = ARKodeSVtolerances(arkode_mem, $(double relTol), tv);
597 if (check_flag(&flag, "ARKodeSVtolerances", 1)) return 1;
598
599 /* Initialize dense matrix data structure and solver */
600 A = SUNDenseMatrix(NEQ, NEQ);
601 if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
602 LS = SUNDenseLinearSolver(y, A);
603 if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;
604
605 /* Attach matrix and linear solver */
606 flag = ARKDlsSetLinearSolver(arkode_mem, LS, A);
607 if (check_flag(&flag, "ARKDlsSetLinearSolver", 1)) return 1;
608
609 /* Set the initial step size if there is one */
610 if ($(int isInitStepSize)) {
611 /* FIXME: We could check if the initial step size is 0 */
612 /* or even NaN and then throw an error */
613 flag = ARKodeSetInitStep(arkode_mem, $(double ss));
614 if (check_flag(&flag, "ARKodeSetInitStep", 1)) return 1;
615 }
616
617 /* Set the Jacobian if there is one */
618 if ($(int isJac)) {
619 flag = ARKDlsSetJacFn(arkode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
620 if (check_flag(&flag, "ARKDlsSetJacFn", 1)) return 1;
621 }
622
623 /* Store initial conditions */
624 for (j = 0; j < NEQ; j++) {
625 ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
626 }
627
628 /* Explicitly set the method */
629 if ($(int method) >= MIN_DIRK_NUM) {
630 flag = ARKodeSetIRKTableNum(arkode_mem, $(int method));
631 if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
632 } else {
633 flag = ARKodeSetERKTableNum(arkode_mem, $(int method));
634 if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
635 }
636
637 /* Main time-stepping loop: calls ARKode to perform the integration */
638 /* Stops when the final time has been reached */
639 for (i = 1; i < $(int nTs); i++) {
640
641 flag = ARKode(arkode_mem, ($vec-ptr:(double *tMut))[i], y, &t, ARK_NORMAL); /* call integrator */
642 if (check_flag(&flag, "ARKode", 1)) break;
643
644 /* Store the results for Haskell */
645 for (j = 0; j < NEQ; j++) {
646 ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
647 }
648
649 /* unsuccessful solve: break */
650 if (flag < 0) {
651 fprintf(stderr,"Solver failure, stopping integration\n");
652 break;
653 }
654 }
655
656 /* Get some final statistics on how the solve progressed */
657
658 flag = ARKodeGetNumSteps(arkode_mem, &nst);
659 check_flag(&flag, "ARKodeGetNumSteps", 1);
660 ($vec-ptr:(long int *diagMut))[0] = nst;
661
662 flag = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
663 check_flag(&flag, "ARKodeGetNumStepAttempts", 1);
664 ($vec-ptr:(long int *diagMut))[1] = nst_a;
665
666 flag = ARKodeGetNumRhsEvals(arkode_mem, &nfe, &nfi);
667 check_flag(&flag, "ARKodeGetNumRhsEvals", 1);
668 ($vec-ptr:(long int *diagMut))[2] = nfe;
669 ($vec-ptr:(long int *diagMut))[3] = nfi;
670
671 flag = ARKodeGetNumLinSolvSetups(arkode_mem, &nsetups);
672 check_flag(&flag, "ARKodeGetNumLinSolvSetups", 1);
673 ($vec-ptr:(long int *diagMut))[4] = nsetups;
674
675 flag = ARKodeGetNumErrTestFails(arkode_mem, &netf);
676 check_flag(&flag, "ARKodeGetNumErrTestFails", 1);
677 ($vec-ptr:(long int *diagMut))[5] = netf;
678
679 flag = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
680 check_flag(&flag, "ARKodeGetNumNonlinSolvIters", 1);
681 ($vec-ptr:(long int *diagMut))[6] = nni;
682
683 flag = ARKodeGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
684 check_flag(&flag, "ARKodeGetNumNonlinSolvConvFails", 1);
685 ($vec-ptr:(long int *diagMut))[7] = ncfn;
686
687 flag = ARKDlsGetNumJacEvals(arkode_mem, &nje);
688 check_flag(&flag, "ARKDlsGetNumJacEvals", 1);
689 ($vec-ptr:(long int *diagMut))[8] = ncfn;
690
691 flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
692 check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
693 ($vec-ptr:(long int *diagMut))[9] = ncfn;
694
695 /* Clean up and return */
696 N_VDestroy(y); /* Free y vector */
697 N_VDestroy(tv); /* Free tv vector */
698 ARKodeFree(&arkode_mem); /* Free integrator memory */
699 SUNLinSolFree(LS); /* Free linear solver */
700 SUNMatDestroy(A); /* Free A matrix */
701
702 return flag;
703 } |]
704 if res == 0
705 then do
706 preD <- V.freeze diagMut
707 let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
708 (fromIntegral $ preD V.!1)
709 (fromIntegral $ preD V.!2)
710 (fromIntegral $ preD V.!3)
711 (fromIntegral $ preD V.!4)
712 (fromIntegral $ preD V.!5)
713 (fromIntegral $ preD V.!6)
714 (fromIntegral $ preD V.!7)
715 (fromIntegral $ preD V.!8)
716 (fromIntegral $ preD V.!9)
717 m <- V.freeze qMatMut
718 return $ Right (m, d)
719 else do
720 return $ Left res
721
722data ButcherTable = ButcherTable { am :: Matrix Double
723 , cv :: Vector Double
724 , bv :: Vector Double
725 , b2v :: Vector Double
726 }
727 deriving Show
728
729data ButcherTable' a = ButcherTable' { am' :: V.Vector a
730 , cv' :: V.Vector a
731 , bv' :: V.Vector a
732 , b2v' :: V.Vector a
733 }
734 deriving Show
735
736butcherTable :: ODEMethod -> ButcherTable
737butcherTable method =
738 case getBT method of
739 Left c -> error $ show c -- FIXME
740 Right (ButcherTable' v w x y, sqp) ->
741 ButcherTable { am = subMatrix (0, 0) (s, s) $ (B.arkSMax >< B.arkSMax) (V.toList v)
742 , cv = subVector 0 s w
743 , bv = subVector 0 s x
744 , b2v = subVector 0 s y
745 }
746 where
747 s = fromIntegral $ sqp V.! 0
748
749getBT :: ODEMethod -> Either Int (ButcherTable' Double, V.Vector Int)
750getBT method = case getButcherTable method of
751 Left c ->
752 Left $ fromIntegral c
753 Right (ButcherTable' a b c d, sqp) ->
754 Right $ ( ButcherTable' (coerce a) (coerce b) (coerce c) (coerce d)
755 , V.map fromIntegral sqp )
756
757getButcherTable :: ODEMethod
758 -> Either CInt (ButcherTable' CDouble, V.Vector CInt)
759getButcherTable method = unsafePerformIO $ do
760 -- ARKode seems to want an ODE in order to set and then get the
761 -- Butcher tableau so here's one to keep it happy
762 let funI :: CDouble -> V.Vector CDouble -> V.Vector CDouble
763 funI _t ys = V.fromList [ ys V.! 0 ]
764 let funE :: CDouble -> V.Vector CDouble -> V.Vector CDouble
765 funE _t ys = V.fromList [ ys V.! 0 ]
766 f0 = V.fromList [ 1.0 ]
767 ts = V.fromList [ 0.0 ]
768 dim = V.length f0
769 nEq :: CLong
770 nEq = fromIntegral dim
771 mN :: CInt
772 mN = fromIntegral $ getMethod method
773
774 btSQP :: V.Vector CInt <- createVector 3
775 btSQPMut <- V.thaw btSQP
776 btAs :: V.Vector CDouble <- createVector (B.arkSMax * B.arkSMax)
777 btAsMut <- V.thaw btAs
778 btCs :: V.Vector CDouble <- createVector B.arkSMax
779 btBs :: V.Vector CDouble <- createVector B.arkSMax
780 btB2s :: V.Vector CDouble <- createVector B.arkSMax
781 btCsMut <- V.thaw btCs
782 btBsMut <- V.thaw btBs
783 btB2sMut <- V.thaw btB2s
784 let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
785 funIOI x y f _ptr = do
786 fImm <- funI x <$> getDataFromContents dim y
787 putDataInContents fImm dim f
788 -- FIXME: I don't understand what this comment means
789 -- Unsafe since the function will be called many times.
790 [CU.exp| int{ 0 } |]
791 let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
792 funIOE x y f _ptr = do
793 fImm <- funE x <$> getDataFromContents dim y
794 putDataInContents fImm dim f
795 -- FIXME: I don't understand what this comment means
796 -- Unsafe since the function will be called many times.
797 [CU.exp| int{ 0 } |]
798 res <- [C.block| int {
799 /* general problem variables */
800
801 int flag; /* reusable error-checking flag */
802 N_Vector y = NULL; /* empty vector for storing solution */
803 void *arkode_mem = NULL; /* empty ARKode memory structure */
804 int i, j; /* reusable loop indices */
805
806 /* general problem parameters */
807
808 realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time */
809 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars */
810
811 /* Initialize data structures */
812
813 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
814 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
815 /* Specify initial condition */
816 for (i = 0; i < NEQ; i++) {
817 NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
818 };
819 arkode_mem = ARKodeCreate(); /* Create the solver memory */
820 if (check_flag((void *)arkode_mem, "ARKodeCreate", 0)) return 1;
821
822 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);
823 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
824
825 if ($(int mN) >= MIN_DIRK_NUM) {
826 flag = ARKodeSetIRKTableNum(arkode_mem, $(int mN));
827 if (check_flag(&flag, "ARKodeSetIRKTableNum", 1)) return 1;
828 } else {
829 flag = ARKodeSetERKTableNum(arkode_mem, $(int mN));
830 if (check_flag(&flag, "ARKodeSetERKTableNum", 1)) return 1;
831 }
832
833 int s, q, p;
834 realtype *ai = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
835 realtype *ae = (realtype *)malloc(ARK_S_MAX * ARK_S_MAX * sizeof(realtype));
836 realtype *ci = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
837 realtype *ce = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
838 realtype *bi = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
839 realtype *be = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
840 realtype *b2i = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
841 realtype *b2e = (realtype *)malloc(ARK_S_MAX * sizeof(realtype));
842 flag = ARKodeGetCurrentButcherTables(arkode_mem, &s, &q, &p, ai, ae, ci, ce, bi, be, b2i, b2e);
843 if (check_flag(&flag, "ARKode", 1)) return 1;
844 $vec-ptr:(int *btSQPMut)[0] = s;
845 $vec-ptr:(int *btSQPMut)[1] = q;
846 $vec-ptr:(int *btSQPMut)[2] = p;
847 for (i = 0; i < s; i++) {
848 for (j = 0; j < s; j++) {
849 /* FIXME: double should be realtype */
850 ($vec-ptr:(double *btAsMut))[i * ARK_S_MAX + j] = ai[i * ARK_S_MAX + j];
851 }
852 }
853
854 for (i = 0; i < s; i++) {
855 ($vec-ptr:(double *btCsMut))[i] = ci[i];
856 ($vec-ptr:(double *btBsMut))[i] = bi[i];
857 ($vec-ptr:(double *btB2sMut))[i] = b2i[i];
858 }
859
860 /* Clean up and return */
861 N_VDestroy(y); /* Free y vector */
862 ARKodeFree(&arkode_mem); /* Free integrator memory */
863
864 return flag;
865 } |]
866 if res == 0
867 then do
868 x <- V.freeze btAsMut
869 y <- V.freeze btSQPMut
870 z <- V.freeze btCsMut
871 u <- V.freeze btBsMut
872 v <- V.freeze btB2sMut
873 return $ Right (ButcherTable' { am' = x, cv' = z, bv' = u, b2v' = v }, y)
874 else do
875 return $ Left res
876
877-- | Adaptive step-size control
878-- functions.
879--
880-- [GSL](https://www.gnu.org/software/gsl/doc/html/ode-initval.html#adaptive-step-size-control)
881-- allows the user to control the step size adjustment using
882-- \(D_i = \epsilon^{abs}s_i + \epsilon^{rel}(a_{y} |y_i| + a_{dy/dt} h |\dot{y}_i|)\) where
883-- \(\epsilon^{abs}\) is the required absolute error, \(\epsilon^{rel}\)
884-- is the required relative error, \(s_i\) is a vector of scaling
885-- factors, \(a_{y}\) is a scaling factor for the solution \(y\) and
886-- \(a_{dydt}\) is a scaling factor for the derivative of the solution \(dy/dt\).
887--
888-- [ARKode](https://computation.llnl.gov/projects/sundials/arkode)
889-- allows the user to control the step size adjustment using
890-- \(\eta^{rel}|y_i| + \eta^{abs}_i\). For compatibility with
891-- [hmatrix-gsl](https://hackage.haskell.org/package/hmatrix-gsl),
892-- tolerances for \(y\) and \(\dot{y}\) can be specified but the latter have no
893-- effect.
894data 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
895 | 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
896 | XX' Double Double Double Double -- ^ include both via relative tolerance
897 -- scaling factors \(a_y\), \(a_{{dy}/{dt}}\); in ARKode terms, the latter is ignored and \(\eta^{rel} = a_{y}\epsilon^{rel}\)
898 | 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 @@
1{-# OPTIONS_GHC -Wall #-}
2
3{-# LANGUAGE QuasiQuotes #-}
4{-# LANGUAGE TemplateHaskell #-}
5{-# LANGUAGE MultiWayIf #-}
6{-# LANGUAGE OverloadedStrings #-}
7{-# LANGUAGE EmptyDataDecls #-}
8
9module Types where
10
11import Foreign.C.Types
12
13import qualified Language.Haskell.TH as TH
14import qualified Language.C.Types as CT
15import qualified Data.Map as Map
16import Language.C.Inline.Context
17
18import qualified Data.Vector.Storable as V
19
20
21data SunVector
22data SunMatrix = SunMatrix { rows :: CInt
23 , cols :: CInt
24 , vals :: V.Vector CDouble
25 }
26
27-- FIXME: Is this true?
28type SunIndexType = CLong
29
30sunTypesTable :: Map.Map CT.TypeSpecifier TH.TypeQ
31sunTypesTable = Map.fromList
32 [
33 (CT.TypeName "sunindextype", [t| SunIndexType |] )
34 , (CT.TypeName "SunVector", [t| SunVector |] )
35 , (CT.TypeName "SunMatrix", [t| SunMatrix |] )
36 ]
37
38sunCtx :: Context
39sunCtx = mempty {ctxTypesTable = sunTypesTable}
40
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 @@
1#include <stdio.h>
2#include <math.h>
3#include <arkode/arkode.h> /* prototypes for ARKODE fcts., consts. */
4#include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
5#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
6#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
7#include <arkode/arkode_direct.h> /* access to ARKDls interface */
8#include <sundials/sundials_types.h> /* definition of type realtype */
9#include <sundials/sundials_math.h>
10
11/* Check function return value...
12 opt == 0 means SUNDIALS function allocates memory so check if
13 returned NULL pointer
14 opt == 1 means SUNDIALS function returns a flag so check if
15 flag >= 0
16 opt == 2 means function allocates memory so check if returned
17 NULL pointer
18*/
19int check_flag(void *flagvalue, const char *funcname, int opt)
20{
21 int *errflag;
22
23 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
24 if (opt == 0 && flagvalue == NULL) {
25 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
26 funcname);
27 return 1; }
28
29 /* Check if flag < 0 */
30 else if (opt == 1) {
31 errflag = (int *) flagvalue;
32 if (*errflag < 0) {
33 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
34 funcname, *errflag);
35 return 1; }}
36
37 /* Check if function returned NULL pointer - no memory allocated */
38 else if (opt == 2 && flagvalue == NULL) {
39 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
40 funcname);
41 return 1; }
42
43 return 0;
44}
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 @@
1/* Check function return value...
2 opt == 0 means SUNDIALS function allocates memory so check if
3 returned NULL pointer
4 opt == 1 means SUNDIALS function returns a flag so check if
5 flag >= 0
6 opt == 2 means function allocates memory so check if returned
7 NULL pointer
8*/
9int 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 @@
1Name: hmatrix-tests 1Name: hmatrix-tests
2Version: 0.6.0.0 2Version: 0.19.0.0
3License: BSD3 3License: BSD3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: 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 @@
1{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} 1{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns -fno-warn-missing-signatures #-}
2{- | 2{- |
3Module : Numeric.GLS.Tests 3Module : Numeric.GLS.Tests
4Copyright : (c) Alberto Ruiz 2014 4Copyright : (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 @@
1{-# LANGUAGE CPP #-} 1{-# LANGUAGE CPP #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} 2{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns -fno-warn-missing-signatures #-}
3{-# LANGUAGE DataKinds #-} 3{-# LANGUAGE DataKinds #-}
4{-# LANGUAGE TypeFamilies #-} 4{-# LANGUAGE TypeFamilies #-}
5{-# LANGUAGE FlexibleContexts #-} 5{-# LANGUAGE FlexibleContexts #-}
@@ -31,7 +31,7 @@ module Numeric.LinearAlgebra.Tests(
31--, runBigTests 31--, runBigTests
32) where 32) where
33 33
34import Numeric.LinearAlgebra hiding (unitary) 34import Numeric.LinearAlgebra
35import Numeric.LinearAlgebra.Devel 35import Numeric.LinearAlgebra.Devel
36import Numeric.LinearAlgebra.Static(L) 36import Numeric.LinearAlgebra.Static(L)
37import Numeric.LinearAlgebra.Tests.Instances 37import Numeric.LinearAlgebra.Tests.Instances
@@ -514,7 +514,7 @@ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
514 514
515-------------------------------------------------------------------------------- 515--------------------------------------------------------------------------------
516 516
517sliceTest = utest "slice test" $ and 517_sliceTest = TestList
518 [ testSlice (chol . trustSym) (gen 5 :: Matrix R) 518 [ testSlice (chol . trustSym) (gen 5 :: Matrix R)
519 , testSlice (chol . trustSym) (gen 5 :: Matrix C) 519 , testSlice (chol . trustSym) (gen 5 :: Matrix C)
520 , testSlice qr (rec :: Matrix R) 520 , testSlice qr (rec :: Matrix R)
@@ -617,7 +617,7 @@ sliceTest = utest "slice test" $ and
617 617
618 test_qrgr n t x = qrgr n (QR x t) 618 test_qrgr n t x = qrgr n (QR x t)
619 619
620 ok_qrgr x = simeq 1E-15 q q' 620 ok_qrgr x = TestCase . assertBool "ok_qrgr" $ simeq 1E-15 q q'
621 where 621 where
622 (q,_) = qr x 622 (q,_) = qr x
623 atau = qrRaw x 623 atau = qrRaw x
@@ -646,7 +646,8 @@ sliceTest = utest "slice test" $ and
646 rec :: Numeric t => Matrix t 646 rec :: Numeric t => Matrix t
647 rec = subMatrix (0,0) (4,5) (gen 5) 647 rec = subMatrix (0,0) (4,5) (gen 5)
648 648
649 testSlice f x@(size->sz@(r,c)) = all (==f x) (map f (g y1 ++ g y2)) 649 testSlice f x@(size->sz@(r,c)) =
650 TestList . map (TestCase . assertEqual "" (f x)) $ (map f (g y1 ++ g y2))
650 where 651 where
651 subm = subMatrix 652 subm = subMatrix
652 g y = [ subm (a*r,b*c) sz y | a <-[0..2], b <- [0..2]] 653 g y = [ subm (a*r,b*c) sz y | a <-[0..2], b <- [0..2]]
@@ -841,7 +842,7 @@ runTests n = do
841 , staticTest 842 , staticTest
842 , intTest 843 , intTest
843 , modularTest 844 , modularTest
844 , sliceTest 845 -- , sliceTest
845 ] 846 ]
846 when (errors c + failures c > 0) exitFailure 847 when (errors c + failures c > 0) exitFailure
847 return () 848 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 @@
1{-# LANGUAGE CPP, FlexibleContexts, UndecidableInstances, FlexibleInstances, ScopedTypeVariables #-} 1{-# LANGUAGE CPP, FlexibleContexts, UndecidableInstances, FlexibleInstances, ScopedTypeVariables #-}
2
3{-# OPTIONS_GHC -fno-warn-orphans #-}
4{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
5
2----------------------------------------------------------------------------- 6-----------------------------------------------------------------------------
3{- | 7{- |
4Module : Numeric.LinearAlgebra.Tests.Instances 8Module : Numeric.LinearAlgebra.Tests.Instances
@@ -62,7 +66,7 @@ instance KnownNat n => Arbitrary (Static.R n) where
62 n :: Int 66 n :: Int
63 n = fromIntegral (natVal (Proxy :: Proxy n)) 67 n = fromIntegral (natVal (Proxy :: Proxy n))
64 68
65 shrink v = [] 69 shrink _v = []
66 70
67instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 71instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
68 arbitrary = do 72 arbitrary = do
@@ -89,7 +93,7 @@ instance (KnownNat n, KnownNat m) => Arbitrary (Static.L m n) where
89 n :: Int 93 n :: Int
90 n = fromIntegral (natVal (Proxy :: Proxy n)) 94 n = fromIntegral (natVal (Proxy :: Proxy n))
91 95
92 shrink mat = [] 96 shrink _mat = []
93 97
94-- a square matrix 98-- a square matrix
95newtype (Sq a) = Sq (Matrix a) deriving Show 99newtype (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 @@
3{-# LANGUAGE TypeFamilies #-} 3{-# LANGUAGE TypeFamilies #-}
4{-# LANGUAGE DataKinds #-} 4{-# LANGUAGE DataKinds #-}
5 5
6{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
7
6----------------------------------------------------------------------------- 8-----------------------------------------------------------------------------
7{- | 9{- |
8Module : Numeric.LinearAlgebra.Tests.Properties 10Module : Numeric.LinearAlgebra.Tests.Properties
@@ -51,14 +53,13 @@ module Numeric.LinearAlgebra.Tests.Properties (
51 , staticVectorBinaryFailProp 53 , staticVectorBinaryFailProp
52) where 54) where
53 55
54import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary) 56import Numeric.LinearAlgebra.HMatrix hiding (Testable)
55import qualified Numeric.LinearAlgebra.Static as Static 57import qualified Numeric.LinearAlgebra.Static as Static
56import Test.QuickCheck 58import Test.QuickCheck
57 59
58import Data.Binary 60import Data.Binary
59import Data.Binary.Get (runGet) 61import Data.Binary.Get (runGet)
60import Data.Either (isLeft) 62import Data.Either (isLeft)
61import Debug.Trace (traceShowId)
62#if MIN_VERSION_base(4,11,0) 63#if MIN_VERSION_base(4,11,0)
63import Prelude hiding ((<>)) 64import Prelude hiding ((<>))
64#endif 65#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 @@
1import Numeric.LinearAlgebra.Tests 1import Numeric.LinearAlgebra.Tests
2 2
3main :: IO ()
3main = runTests 20 4main = 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 @@
1import Numeric.GSL.Tests 1import Numeric.GSL.Tests
2 2
3main :: IO ()
3main = runTests 20 4main = runTests 20