diff options
author | Peter Dobsan <pdobsan@gmail.com> | 2018-05-03 20:08:33 +0200 |
---|---|---|
committer | Peter Dobsan <pdobsan@gmail.com> | 2018-05-03 20:08:33 +0200 |
commit | cafdc664c01ea7392c81c352b5c5444dc2963531 (patch) | |
tree | c6d9a758fa7c36730c0468b393a6dc8c47cbfac2 /packages/sundials/src/Numeric/Sundials/Arkode.hsc | |
parent | ea1bfea4486f8f2c646f82dabd1ff9a222b68506 (diff) | |
parent | 1675813d8f540af9832a78c7a7a40bbdf1cec42c (diff) |
Merge remote-tracking branch 'upstream/master'
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/Arkode.hsc')
-rw-r--r-- | packages/sundials/src/Numeric/Sundials/Arkode.hsc | 204 |
1 files changed, 204 insertions, 0 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/Arkode.hsc b/packages/sundials/src/Numeric/Sundials/Arkode.hsc new file mode 100644 index 0000000..0850258 --- /dev/null +++ b/packages/sundials/src/Numeric/Sundials/Arkode.hsc | |||
@@ -0,0 +1,204 @@ | |||
1 | {-# LANGUAGE QuasiQuotes #-} | ||
2 | {-# LANGUAGE TemplateHaskell #-} | ||
3 | {-# LANGUAGE OverloadedStrings #-} | ||
4 | {-# LANGUAGE EmptyDataDecls #-} | ||
5 | |||
6 | module Numeric.Sundials.Arkode where | ||
7 | |||
8 | import Foreign | ||
9 | import Foreign.C.Types | ||
10 | |||
11 | import Language.C.Types as CT | ||
12 | |||
13 | import qualified Data.Vector.Storable as VS | ||
14 | import qualified Data.Vector.Storable.Mutable as VM | ||
15 | |||
16 | import qualified Language.Haskell.TH as TH | ||
17 | import qualified Data.Map as Map | ||
18 | import Language.C.Inline.Context | ||
19 | |||
20 | import qualified Data.Vector.Storable as V | ||
21 | |||
22 | |||
23 | #include <stdio.h> | ||
24 | #include <sundials/sundials_nvector.h> | ||
25 | #include <sundials/sundials_matrix.h> | ||
26 | #include <nvector/nvector_serial.h> | ||
27 | #include <sunmatrix/sunmatrix_dense.h> | ||
28 | #include <arkode/arkode.h> | ||
29 | #include <cvode/cvode.h> | ||
30 | |||
31 | |||
32 | data SunVector | ||
33 | data SunMatrix = SunMatrix { rows :: CInt | ||
34 | , cols :: CInt | ||
35 | , vals :: V.Vector CDouble | ||
36 | } | ||
37 | |||
38 | -- | This is true only if configured/ built as 64 bits | ||
39 | type SunIndexType = CLong | ||
40 | |||
41 | sunTypesTable :: Map.Map TypeSpecifier TH.TypeQ | ||
42 | sunTypesTable = Map.fromList | ||
43 | [ | ||
44 | (TypeName "sunindextype", [t| SunIndexType |] ) | ||
45 | , (TypeName "SunVector", [t| SunVector |] ) | ||
46 | , (TypeName "SunMatrix", [t| SunMatrix |] ) | ||
47 | ] | ||
48 | |||
49 | sunCtx :: Context | ||
50 | sunCtx = mempty {ctxTypesTable = sunTypesTable} | ||
51 | |||
52 | getMatrixDataFromContents :: Ptr SunMatrix -> IO SunMatrix | ||
53 | getMatrixDataFromContents ptr = do | ||
54 | qtr <- getContentMatrixPtr ptr | ||
55 | rs <- getNRows qtr | ||
56 | cs <- getNCols qtr | ||
57 | rtr <- getMatrixData qtr | ||
58 | vs <- vectorFromC (fromIntegral $ rs * cs) rtr | ||
59 | return $ SunMatrix { rows = rs, cols = cs, vals = vs } | ||
60 | |||
61 | putMatrixDataFromContents :: SunMatrix -> Ptr SunMatrix -> IO () | ||
62 | putMatrixDataFromContents mat ptr = do | ||
63 | let rs = rows mat | ||
64 | cs = cols mat | ||
65 | vs = vals mat | ||
66 | qtr <- getContentMatrixPtr ptr | ||
67 | putNRows rs qtr | ||
68 | putNCols cs qtr | ||
69 | rtr <- getMatrixData qtr | ||
70 | vectorToC vs (fromIntegral $ rs * cs) rtr | ||
71 | |||
72 | instance Storable SunMatrix where | ||
73 | poke = flip putMatrixDataFromContents | ||
74 | peek = getMatrixDataFromContents | ||
75 | sizeOf _ = error "sizeOf not supported for SunMatrix" | ||
76 | alignment _ = error "alignment not supported for SunMatrix" | ||
77 | |||
78 | vectorFromC :: Storable a => Int -> Ptr a -> IO (VS.Vector a) | ||
79 | vectorFromC len ptr = do | ||
80 | ptr' <- newForeignPtr_ ptr | ||
81 | VS.freeze $ VM.unsafeFromForeignPtr0 ptr' len | ||
82 | |||
83 | vectorToC :: Storable a => VS.Vector a -> Int -> Ptr a -> IO () | ||
84 | vectorToC vec len ptr = do | ||
85 | ptr' <- newForeignPtr_ ptr | ||
86 | VS.copy (VM.unsafeFromForeignPtr0 ptr' len) vec | ||
87 | |||
88 | getDataFromContents :: Int -> Ptr SunVector -> IO (VS.Vector CDouble) | ||
89 | getDataFromContents len ptr = do | ||
90 | qtr <- getContentPtr ptr | ||
91 | rtr <- getData qtr | ||
92 | vectorFromC len rtr | ||
93 | |||
94 | putDataInContents :: Storable a => VS.Vector a -> Int -> Ptr b -> IO () | ||
95 | putDataInContents vec len ptr = do | ||
96 | qtr <- getContentPtr ptr | ||
97 | rtr <- getData qtr | ||
98 | vectorToC vec len rtr | ||
99 | |||
100 | #def typedef struct _generic_N_Vector SunVector; | ||
101 | #def typedef struct _N_VectorContent_Serial SunContent; | ||
102 | |||
103 | #def typedef struct _generic_SUNMatrix SunMatrix; | ||
104 | #def typedef struct _SUNMatrixContent_Dense SunMatrixContent; | ||
105 | |||
106 | getContentMatrixPtr :: Storable a => Ptr b -> IO a | ||
107 | getContentMatrixPtr ptr = (#peek SunMatrix, content) ptr | ||
108 | |||
109 | getNRows :: Ptr b -> IO CInt | ||
110 | getNRows ptr = (#peek SunMatrixContent, M) ptr | ||
111 | putNRows :: CInt -> Ptr b -> IO () | ||
112 | putNRows nr ptr = (#poke SunMatrixContent, M) ptr nr | ||
113 | |||
114 | getNCols :: Ptr b -> IO CInt | ||
115 | getNCols ptr = (#peek SunMatrixContent, N) ptr | ||
116 | putNCols :: CInt -> Ptr b -> IO () | ||
117 | putNCols nc ptr = (#poke SunMatrixContent, N) ptr nc | ||
118 | |||
119 | getMatrixData :: Storable a => Ptr b -> IO a | ||
120 | getMatrixData ptr = (#peek SunMatrixContent, data) ptr | ||
121 | |||
122 | getContentPtr :: Storable a => Ptr b -> IO a | ||
123 | getContentPtr ptr = (#peek SunVector, content) ptr | ||
124 | |||
125 | getData :: Storable a => Ptr b -> IO a | ||
126 | getData ptr = (#peek SunContent, data) ptr | ||
127 | |||
128 | cV_ADAMS :: Int | ||
129 | cV_ADAMS = #const CV_ADAMS | ||
130 | cV_BDF :: Int | ||
131 | cV_BDF = #const CV_BDF | ||
132 | |||
133 | arkSMax :: Int | ||
134 | arkSMax = #const ARK_S_MAX | ||
135 | |||
136 | mIN_DIRK_NUM, mAX_DIRK_NUM :: Int | ||
137 | mIN_DIRK_NUM = #const MIN_DIRK_NUM | ||
138 | mAX_DIRK_NUM = #const MAX_DIRK_NUM | ||
139 | |||
140 | -- FIXME: We could just use inline-c instead | ||
141 | |||
142 | -- Butcher table accessors -- implicit | ||
143 | sDIRK_2_1_2 :: Int | ||
144 | sDIRK_2_1_2 = #const SDIRK_2_1_2 | ||
145 | bILLINGTON_3_3_2 :: Int | ||
146 | bILLINGTON_3_3_2 = #const BILLINGTON_3_3_2 | ||
147 | tRBDF2_3_3_2 :: Int | ||
148 | tRBDF2_3_3_2 = #const TRBDF2_3_3_2 | ||
149 | kVAERNO_4_2_3 :: Int | ||
150 | kVAERNO_4_2_3 = #const KVAERNO_4_2_3 | ||
151 | aRK324L2SA_DIRK_4_2_3 :: Int | ||
152 | aRK324L2SA_DIRK_4_2_3 = #const ARK324L2SA_DIRK_4_2_3 | ||
153 | cASH_5_2_4 :: Int | ||
154 | cASH_5_2_4 = #const CASH_5_2_4 | ||
155 | cASH_5_3_4 :: Int | ||
156 | cASH_5_3_4 = #const CASH_5_3_4 | ||
157 | sDIRK_5_3_4 :: Int | ||
158 | sDIRK_5_3_4 = #const SDIRK_5_3_4 | ||
159 | kVAERNO_5_3_4 :: Int | ||
160 | kVAERNO_5_3_4 = #const KVAERNO_5_3_4 | ||
161 | aRK436L2SA_DIRK_6_3_4 :: Int | ||
162 | aRK436L2SA_DIRK_6_3_4 = #const ARK436L2SA_DIRK_6_3_4 | ||
163 | kVAERNO_7_4_5 :: Int | ||
164 | kVAERNO_7_4_5 = #const KVAERNO_7_4_5 | ||
165 | aRK548L2SA_DIRK_8_4_5 :: Int | ||
166 | aRK548L2SA_DIRK_8_4_5 = #const ARK548L2SA_DIRK_8_4_5 | ||
167 | |||
168 | -- #define DEFAULT_DIRK_2 SDIRK_2_1_2 | ||
169 | -- #define DEFAULT_DIRK_3 ARK324L2SA_DIRK_4_2_3 | ||
170 | -- #define DEFAULT_DIRK_4 SDIRK_5_3_4 | ||
171 | -- #define DEFAULT_DIRK_5 ARK548L2SA_DIRK_8_4_5 | ||
172 | |||
173 | -- Butcher table accessors -- explicit | ||
174 | hEUN_EULER_2_1_2 :: Int | ||
175 | hEUN_EULER_2_1_2 = #const HEUN_EULER_2_1_2 | ||
176 | bOGACKI_SHAMPINE_4_2_3 :: Int | ||
177 | bOGACKI_SHAMPINE_4_2_3 = #const BOGACKI_SHAMPINE_4_2_3 | ||
178 | aRK324L2SA_ERK_4_2_3 :: Int | ||
179 | aRK324L2SA_ERK_4_2_3 = #const ARK324L2SA_ERK_4_2_3 | ||
180 | zONNEVELD_5_3_4 :: Int | ||
181 | zONNEVELD_5_3_4 = #const ZONNEVELD_5_3_4 | ||
182 | aRK436L2SA_ERK_6_3_4 :: Int | ||
183 | aRK436L2SA_ERK_6_3_4 = #const ARK436L2SA_ERK_6_3_4 | ||
184 | sAYFY_ABURUB_6_3_4 :: Int | ||
185 | sAYFY_ABURUB_6_3_4 = #const SAYFY_ABURUB_6_3_4 | ||
186 | cASH_KARP_6_4_5 :: Int | ||
187 | cASH_KARP_6_4_5 = #const CASH_KARP_6_4_5 | ||
188 | fEHLBERG_6_4_5 :: Int | ||
189 | fEHLBERG_6_4_5 = #const FEHLBERG_6_4_5 | ||
190 | dORMAND_PRINCE_7_4_5 :: Int | ||
191 | dORMAND_PRINCE_7_4_5 = #const DORMAND_PRINCE_7_4_5 | ||
192 | aRK548L2SA_ERK_8_4_5 :: Int | ||
193 | aRK548L2SA_ERK_8_4_5 = #const ARK548L2SA_ERK_8_4_5 | ||
194 | vERNER_8_5_6 :: Int | ||
195 | vERNER_8_5_6 = #const VERNER_8_5_6 | ||
196 | fEHLBERG_13_7_8 :: Int | ||
197 | fEHLBERG_13_7_8 = #const FEHLBERG_13_7_8 | ||
198 | |||
199 | -- #define DEFAULT_ERK_2 HEUN_EULER_2_1_2 | ||
200 | -- #define DEFAULT_ERK_3 BOGACKI_SHAMPINE_4_2_3 | ||
201 | -- #define DEFAULT_ERK_4 ZONNEVELD_5_3_4 | ||
202 | -- #define DEFAULT_ERK_5 CASH_KARP_6_4_5 | ||
203 | -- #define DEFAULT_ERK_6 VERNER_8_5_6 | ||
204 | -- #define DEFAULT_ERK_8 FEHLBERG_13_7_8 | ||