summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/CVode/ODE.hs')
-rw-r--r--packages/sundials/src/Numeric/Sundials/CVode/ODE.hs41
1 files changed, 21 insertions, 20 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
index 1cd072f..159fbe2 100644
--- a/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/CVode/ODE.hs
@@ -79,6 +79,7 @@ import Data.Maybe (isJust)
79 79
80import Foreign.C.Types (CDouble, CInt, CLong) 80import Foreign.C.Types (CDouble, CInt, CLong)
81import Foreign.Ptr (Ptr) 81import Foreign.Ptr (Ptr)
82import Foreign.Storable (poke)
82 83
83import qualified Data.Vector.Storable as V 84import qualified Data.Vector.Storable as V
84 85
@@ -90,10 +91,10 @@ import Numeric.LinearAlgebra.Devel (createVector)
90import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows, 91import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
91 cols, toLists, size, reshape) 92 cols, toLists, size, reshape)
92 93
93import qualified Numeric.Sundials.CLangToHaskellTypes as T 94import Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF,
94import Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF) 95 getDataFromContents, putDataInContents)
95import Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian) 96import qualified Numeric.Sundials.Arkode as T
96import qualified Numeric.Sundials.ODEOpts as SO 97import Numeric.Sundials.ODEOpts (ODEOpts(..), Jacobian, SundialsDiagnostics(..))
97 98
98 99
99C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx) 100C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
@@ -195,7 +196,7 @@ odeSolveVWith' ::
195 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 196 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
196 -> V.Vector Double -- ^ Initial conditions 197 -> V.Vector Double -- ^ Initial conditions
197 -> V.Vector Double -- ^ Desired solution times 198 -> V.Vector Double -- ^ Desired solution times
198 -> Either Int (Matrix Double, SO.SundialsDiagnostics) -- ^ Error code or solution 199 -> Either Int (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
199odeSolveVWith' opts method control initStepSize f y0 tt = 200odeSolveVWith' opts method control initStepSize f y0 tt =
200 case solveOdeC (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts) 201 case solveOdeC (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
201 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control) 202 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
@@ -229,7 +230,7 @@ solveOdeC ::
229 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 230 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
230 -> V.Vector CDouble -- ^ Initial conditions 231 -> V.Vector CDouble -- ^ Initial conditions
231 -> V.Vector CDouble -- ^ Desired solution times 232 -> V.Vector CDouble -- ^ Desired solution times
232 -> Either CInt ((V.Vector CDouble), SO.SundialsDiagnostics) -- ^ Error code or solution 233 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
233solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts = 234solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts =
234 unsafePerformIO $ do 235 unsafePerformIO $ do
235 236
@@ -257,9 +258,9 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
257 funIO x y f _ptr = do 258 funIO x y f _ptr = do
258 -- Convert the pointer we get from C (y) to a vector, and then 259 -- Convert the pointer we get from C (y) to a vector, and then
259 -- apply the user-supplied function. 260 -- apply the user-supplied function.
260 fImm <- fun x <$> SO.getDataFromContents dim y 261 fImm <- fun x <$> getDataFromContents dim y
261 -- Fill in the provided pointer with the resulting vector. 262 -- Fill in the provided pointer with the resulting vector.
262 SO.putDataInContents fImm dim f 263 putDataInContents fImm dim f
263 -- FIXME: I don't understand what this comment means 264 -- FIXME: I don't understand what this comment means
264 -- Unsafe since the function will be called many times. 265 -- Unsafe since the function will be called many times.
265 [CU.exp| int{ 0 } |] 266 [CU.exp| int{ 0 } |]
@@ -271,8 +272,8 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
271 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do 272 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
272 case jacH of 273 case jacH of
273 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined" 274 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
274 Just jacI -> do j <- jacI t <$> SO.getDataFromContents dim y 275 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
275 SO.putMatrixDataFromContents j jacS 276 poke jacS j
276 -- FIXME: I don't understand what this comment means 277 -- FIXME: I don't understand what this comment means
277 -- Unsafe since the function will be called many times. 278 -- Unsafe since the function will be called many times.
278 [CU.exp| int{ 0 } |] 279 [CU.exp| int{ 0 } |]
@@ -431,16 +432,16 @@ solveOdeC maxNumSteps_ minStep_ method initStepSize jacH (aTols, rTol) fun f0 ts
431 if res == 0 432 if res == 0
432 then do 433 then do
433 preD <- V.freeze diagMut 434 preD <- V.freeze diagMut
434 let d = SO.SundialsDiagnostics (fromIntegral $ preD V.!0) 435 let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
435 (fromIntegral $ preD V.!1) 436 (fromIntegral $ preD V.!1)
436 (fromIntegral $ preD V.!2) 437 (fromIntegral $ preD V.!2)
437 (fromIntegral $ preD V.!3) 438 (fromIntegral $ preD V.!3)
438 (fromIntegral $ preD V.!4) 439 (fromIntegral $ preD V.!4)
439 (fromIntegral $ preD V.!5) 440 (fromIntegral $ preD V.!5)
440 (fromIntegral $ preD V.!6) 441 (fromIntegral $ preD V.!6)
441 (fromIntegral $ preD V.!7) 442 (fromIntegral $ preD V.!7)
442 (fromIntegral $ preD V.!8) 443 (fromIntegral $ preD V.!8)
443 (fromIntegral $ preD V.!9) 444 (fromIntegral $ preD V.!9)
444 m <- V.freeze qMatMut 445 m <- V.freeze qMatMut
445 return $ Right (m, d) 446 return $ Right (m, d)
446 else do 447 else do