summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs')
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs56
1 files changed, 40 insertions, 16 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index a8d418b..13b7eb8 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -125,6 +125,7 @@ import Data.Maybe (isJust)
125 125
126import Foreign.C.Types (CDouble, CInt, CLong) 126import Foreign.C.Types (CDouble, CInt, CLong)
127import Foreign.Ptr (Ptr) 127import Foreign.Ptr (Ptr)
128import Foreign.Storable (poke)
128 129
129import qualified Data.Vector.Storable as V 130import qualified Data.Vector.Storable as V
130 131
@@ -139,10 +140,33 @@ import Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, (><),
139 subMatrix, rows, cols, toLists, 140 subMatrix, rows, cols, toLists,
140 size, subVector) 141 size, subVector)
141 142
142import qualified Numeric.Sundials.CLangToHaskellTypes as T
143import Numeric.Sundials.Arkode
144import qualified Numeric.Sundials.Arkode as B
145import qualified Numeric.Sundials.ODEOpts as SO 143import qualified Numeric.Sundials.ODEOpts as SO
144import qualified Numeric.Sundials.Arkode as T
145import Numeric.Sundials.Arkode (getDataFromContents, putDataInContents, arkSMax,
146 sDIRK_2_1_2,
147 bILLINGTON_3_3_2,
148 tRBDF2_3_3_2,
149 kVAERNO_4_2_3,
150 aRK324L2SA_DIRK_4_2_3,
151 cASH_5_2_4,
152 cASH_5_3_4,
153 sDIRK_5_3_4,
154 kVAERNO_5_3_4,
155 aRK436L2SA_DIRK_6_3_4,
156 kVAERNO_7_4_5,
157 aRK548L2SA_DIRK_8_4_5,
158 hEUN_EULER_2_1_2,
159 bOGACKI_SHAMPINE_4_2_3,
160 aRK324L2SA_ERK_4_2_3,
161 zONNEVELD_5_3_4,
162 aRK436L2SA_ERK_6_3_4,
163 sAYFY_ABURUB_6_3_4,
164 cASH_KARP_6_4_5,
165 fEHLBERG_6_4_5,
166 dORMAND_PRINCE_7_4_5,
167 aRK548L2SA_ERK_8_4_5,
168 vERNER_8_5_6,
169 fEHLBERG_13_7_8)
146 170
147 171
148C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx) 172C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)
@@ -451,9 +475,9 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
451 funIO x y f _ptr = do 475 funIO x y f _ptr = do
452 -- Convert the pointer we get from C (y) to a vector, and then 476 -- Convert the pointer we get from C (y) to a vector, and then
453 -- apply the user-supplied function. 477 -- apply the user-supplied function.
454 fImm <- fun x <$> SO.getDataFromContents dim y 478 fImm <- fun x <$> getDataFromContents dim y
455 -- Fill in the provided pointer with the resulting vector. 479 -- Fill in the provided pointer with the resulting vector.
456 SO.putDataInContents fImm dim f 480 putDataInContents fImm dim f
457 -- FIXME: I don't understand what this comment means 481 -- FIXME: I don't understand what this comment means
458 -- Unsafe since the function will be called many times. 482 -- Unsafe since the function will be called many times.
459 [CU.exp| int{ 0 } |] 483 [CU.exp| int{ 0 } |]
@@ -465,8 +489,8 @@ solveOdeC method initStepSize jacH (absTols, relTol) fun f0 ts = unsafePerformIO
465 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do 489 jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
466 case jacH of 490 case jacH of
467 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined" 491 Nothing -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
468 Just jacI -> do j <- jacI t <$> SO.getDataFromContents dim y 492 Just jacI -> do j <- jacI t <$> getDataFromContents dim y
469 SO.putMatrixDataFromContents j jacS 493 poke jacS j
470 -- FIXME: I don't understand what this comment means 494 -- FIXME: I don't understand what this comment means
471 -- Unsafe since the function will be called many times. 495 -- Unsafe since the function will be called many times.
472 [CU.exp| int{ 0 } |] 496 [CU.exp| int{ 0 } |]
@@ -675,7 +699,7 @@ butcherTable method =
675 case getBT method of 699 case getBT method of
676 Left c -> error $ show c -- FIXME 700 Left c -> error $ show c -- FIXME
677 Right (ButcherTable' v w x y, sqp) -> 701 Right (ButcherTable' v w x y, sqp) ->
678 ButcherTable { am = subMatrix (0, 0) (s, s) $ (B.arkSMax >< B.arkSMax) (V.toList v) 702 ButcherTable { am = subMatrix (0, 0) (s, s) $ (arkSMax >< arkSMax) (V.toList v)
679 , cv = subVector 0 s w 703 , cv = subVector 0 s w
680 , bv = subVector 0 s x 704 , bv = subVector 0 s x
681 , b2v = subVector 0 s y 705 , b2v = subVector 0 s y
@@ -710,25 +734,25 @@ getButcherTable method = unsafePerformIO $ do
710 734
711 btSQP :: V.Vector CInt <- createVector 3 735 btSQP :: V.Vector CInt <- createVector 3
712 btSQPMut <- V.thaw btSQP 736 btSQPMut <- V.thaw btSQP
713 btAs :: V.Vector CDouble <- createVector (B.arkSMax * B.arkSMax) 737 btAs :: V.Vector CDouble <- createVector (arkSMax * arkSMax)
714 btAsMut <- V.thaw btAs 738 btAsMut <- V.thaw btAs
715 btCs :: V.Vector CDouble <- createVector B.arkSMax 739 btCs :: V.Vector CDouble <- createVector arkSMax
716 btBs :: V.Vector CDouble <- createVector B.arkSMax 740 btBs :: V.Vector CDouble <- createVector arkSMax
717 btB2s :: V.Vector CDouble <- createVector B.arkSMax 741 btB2s :: V.Vector CDouble <- createVector arkSMax
718 btCsMut <- V.thaw btCs 742 btCsMut <- V.thaw btCs
719 btBsMut <- V.thaw btBs 743 btBsMut <- V.thaw btBs
720 btB2sMut <- V.thaw btB2s 744 btB2sMut <- V.thaw btB2s
721 let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt 745 let funIOI :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
722 funIOI x y f _ptr = do 746 funIOI x y f _ptr = do
723 fImm <- funI x <$> SO.getDataFromContents dim y 747 fImm <- funI x <$> getDataFromContents dim y
724 SO.putDataInContents fImm dim f 748 putDataInContents fImm dim f
725 -- FIXME: I don't understand what this comment means 749 -- FIXME: I don't understand what this comment means
726 -- Unsafe since the function will be called many times. 750 -- Unsafe since the function will be called many times.
727 [CU.exp| int{ 0 } |] 751 [CU.exp| int{ 0 } |]
728 let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt 752 let funIOE :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
729 funIOE x y f _ptr = do 753 funIOE x y f _ptr = do
730 fImm <- funE x <$> SO.getDataFromContents dim y 754 fImm <- funE x <$> getDataFromContents dim y
731 SO.putDataInContents fImm dim f 755 putDataInContents fImm dim f
732 -- FIXME: I don't understand what this comment means 756 -- FIXME: I don't understand what this comment means
733 -- Unsafe since the function will be called many times. 757 -- Unsafe since the function will be called many times.
734 [CU.exp| int{ 0 } |] 758 [CU.exp| int{ 0 } |]