summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL/Root.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL/Root.hs')
-rw-r--r--lib/Numeric/GSL/Root.hs44
1 files changed, 37 insertions, 7 deletions
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index cd2982a..25ec5f5 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -45,6 +45,7 @@ main = do
45----------------------------------------------------------------------------- 45-----------------------------------------------------------------------------
46 46
47module Numeric.GSL.Root ( 47module Numeric.GSL.Root (
48 uniRoot, UniRootMethod(..),
48 root, RootMethod(..), 49 root, RootMethod(..),
49 rootJ, RootMethodJ(..), 50 rootJ, RootMethodJ(..),
50) where 51) where
@@ -58,6 +59,36 @@ import System.IO.Unsafe(unsafePerformIO)
58 59
59------------------------------------------------------------------------- 60-------------------------------------------------------------------------
60 61
62data UniRootMethod = Bisection
63 | FalsePos
64 | Brent
65 deriving (Enum, Eq, Show, Bounded)
66
67uniRoot :: UniRootMethod
68 -> Double
69 -> Int
70 -> (Double -> Double)
71 -> Double
72 -> Double
73 -> (Double, Matrix Double)
74uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
75
76uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
77 fp <- mkDoublefun f
78 rawpath <- createMIO maxit 4
79 (c_root m fp epsrel (fi maxit) xl xu)
80 "root"
81 let it = round (rawpath @@> (maxit-1,0))
82 path = takeRows it rawpath
83 [sol] = toLists $ dropRows (it-1) path
84 freeHaskellFunPtr fp
85 return (sol !! 1, path)
86
87foreign import ccall safe "root"
88 c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM
89
90-------------------------------------------------------------------------
91
61data RootMethod = Hybrids 92data RootMethod = Hybrids
62 | Hybrid 93 | Hybrid
63 | DNewton 94 | DNewton
@@ -82,8 +113,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
82 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 113 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
83 rawpath <- vec xiv $ \xiv' -> 114 rawpath <- vec xiv $ \xiv' ->
84 createMIO maxit (2*n+1) 115 createMIO maxit (2*n+1)
85 (c_root m fp epsabs (fi maxit) // xiv') 116 (c_multiroot m fp epsabs (fi maxit) // xiv')
86 "root" 117 "multiroot"
87 let it = round (rawpath @@> (maxit-1,0)) 118 let it = round (rawpath @@> (maxit-1,0))
88 path = takeRows it rawpath 119 path = takeRows it rawpath
89 [sol] = toLists $ dropRows (it-1) path 120 [sol] = toLists $ dropRows (it-1) path
@@ -91,8 +122,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
91 return (take n $ drop 1 sol, path) 122 return (take n $ drop 1 sol, path)
92 123
93 124
94foreign import ccall safe "root" 125foreign import ccall safe "multiroot"
95 c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM 126 c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
96 127
97------------------------------------------------------------------------- 128-------------------------------------------------------------------------
98 129
@@ -121,7 +152,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
121 rawpath <- vec xiv $ \xiv' -> 152 rawpath <- vec xiv $ \xiv' ->
122 createMIO maxit (2*n+1) 153 createMIO maxit (2*n+1)
123 (c_rootj m fp jp epsabs (fi maxit) // xiv') 154 (c_rootj m fp jp epsabs (fi maxit) // xiv')
124 "root" 155 "multiroot"
125 let it = round (rawpath @@> (maxit-1,0)) 156 let it = round (rawpath @@> (maxit-1,0))
126 path = takeRows it rawpath 157 path = takeRows it rawpath
127 [sol] = toLists $ dropRows (it-1) path 158 [sol] = toLists $ dropRows (it-1) path
@@ -129,8 +160,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
129 freeHaskellFunPtr jp 160 freeHaskellFunPtr jp
130 return (take n $ drop 1 sol, path) 161 return (take n $ drop 1 sol, path)
131 162
132 163foreign import ccall safe "multirootj"
133foreign import ccall safe "rootj"
134 c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM 164 c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
135 165
136------------------------------------------------------- 166-------------------------------------------------------